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1.0  INTRODUCTION 

During  the  past  few  years,  improvements  made  on  the  AEDC  von  Kdrmdn  Gas 
Dynamics  Facility  (VKF)  three-degree-of-freedom  (3DOF)  test  mechanism  have  led  to 
greater  usage  and  increased  emphasis  on  data  quality  and  data  reduction  techniques.  These 
improvements  consist  primarily  of  the  addition  to  the  original  system  as  reported  in  Ref, 

1 of  roll  position  measurement  and  increased  axial  and  normal  load  capacity. 

Before  the  system  alterations  were  made,  only  a limited  amount  of  data  was  reduced 
to  final  form.  This  was  due  primarily  to  the  relatively  large  computer  time  involved  in 
the  analysis  of  data  and  the  limitations  imposed  by  the  data  reduction  procedure.  The 
analysis  of  data  for  the  extraction  of  meaningful  aerodynamic  coefficients  was  hampered 
by  the  use  of  data  reduction  procedures  based  on  linear  aerodynamics  and  on  mass  and 
aerodynamic  symmetry  properties,  of  the  model  (Ref.  2).  Since  the  majority  of  models 
tested  fail  to  conform  to  either  one  or  both  of  these  restrictions,  more  exact  data  reduction 
techniques  are  required. 

The  capability  of  the  3DOF  system  to  measure  the  roU  position  accurately  permits 
data  to  be  analyzed  in  the  body-axis  system.  Use  of  this  system  allows  analysis  of  the 
effects  of  model  inertia  anomalies  and  asymmetric  aerodynamics  on  3 DOF  motion. 
Moreover,  increased  emphasis  has  been  placed  on  studying  the  effects  of  nonlinear 
aerodynamics  for  both  symmetric  (Ref.  3)  and  asymmetric  configurations  (Ref.  4). 

This  report  presents  improved  data  reduction  and  analysis  procedures  so  that  the 
increased  capabilities  of  the  3 DOF  system  may  be  utilized.  A modified  Newton-Raphson 
technique  (hereafter  referred  to  as  the  Chapman-Kirk  method,  Ref.  5)  was  employed  to 
develop  two  programs  designed  to  extract  aerodynamic  coefficients  from  3 DOF  data.  The 
asymmetric  and  symmetric  programs  utilize  the  full  dynamic  and  kinematic  equations  and 
assume  mirror  plane  symmetry  and  axial  symmetry,  respectively.  Other  data  reduction 
methods  currently  used  in  aerodynamic  parameter  extraction  are  discussed  in  Section  2.0. 
The  formulation  of  the  mathematical  aspects  of  the  Chapman-Kirk  method  and  the 
development  of  the  governing  3 DOF  equations  to  be  used  in  the  extraction  methods  are 
presented  in  Section  3.0.  Section  4.0  gives  a detailed  description  of  the  resulting  programs. 
In  Section  5.0,  the  validity  of  the  improved  coefficient  extraction  technique  is  studied 
using  computer-generated  motion  and  3DOF  data  from  Tunnel  A.  In  addition,  bench  test 
data  are  used  to  support  error  estimates  of  the  extracted  coefficients. 

2.0  METHODS  OF  3DOF  ANALYSIS 

During  the  past  20  years,  much  progress  has  been  made  in  parameter  identification 
techniques.  Early  methods  were  difficult  to  implement  and  time  consuming  in  application. 
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The  need  for  improved  techniques  was  shown  by  the  disagreement  found  between 
extrapolated  wind  tunnel  and  flight  test  results  (Ref.  6).  The  following  discussion 
concentrates  on  several  types  of  identification  techniques  currently  used  in  processing  data, 
both  in  the  wind  tunnel  and  in  full  scale.  The  advantages  and  drawbacks  of  each  type 
are  addressed. 

Parameter  extraction  methods  may  be  broadly  classifled  into  three  categories:  1)  the 
so-called  equation  error  methods,  2)  output  error  methods,  and  3)  advanced  methods. 
The  equation  error  methods  assume  a performance  criterion  (i.e.,  a set  of  differential  or 
algebraic  constraints  on  the  vehicle  aerodynamics  and/or  motion)  that  minimizes  the  square 
of  the  equation  error.  All  methods  that  fall  into  this  group  are  basically  least  squares 
methods,  and  generally  all  response  variables  and  derivatives  must  be  measured.  Application 
of  this  method  gives  a set  of  n linear  equations  with  n unknowns.  Since  the  equations 
to  be  solved  are  linear,  the  equation  error  methods  are  single-step  processes.  As  such, 
they  are  normally  used  only  as  a startup  method  for  output  error  methods.  In  addition, 
equation  error  methods  are  known  to  give  biased  estimates  in  the  presence  of  measurement 
noise. 

Output  error  methods  minimize  the  square  of  the  error  between  the  actual  system 
output  and  the  modeled  output.  Characteristic  of  these  methods  is  the  ability  to  extract 
meaningful  parameters  in  the  presence  of  measurement  noise.  In  addition,  process  noise 
(modeling  errors)  produces  biased  results  in  the  extracted  coefficients  and  may  result  in 
lack  of  convergence  of  the  iterative  scheme.  Typical  output  error  methods  are 
Newton-Raphson,  gradient,  Kalman  filter,  and  modified  Newton-Raphson.  The  modified 
Newton-Raphson  method,  also  referred  to  as  the  differential  correction  or  quasilinearization 
method,  can  be  shown  to  be  identical  to  the  Kalman  filter  method  for  the  case  of  no 
measurement  noise. 

The  maximum  likelihood  method  (as  described  in  Ref.  7)  is  the  most  advanced 
parameter  identification  technique  currently  available  and  is  capable  of  handling  both 
process  and  measurement  noise.  This  technique  may  be  subdivided  into  three  basic  steps. 
The  first  step  employs  the  Kalman  filter  to  estimate  the  magnitude  of  the  states  and 
to  generate  a residual  sequence.  A modified  Newton-Raphson  algorithm  is  then  used  for 
the  parameter  estimates.  Tire  final  step  uses  an  algorithm  to  estimate  the  noise  statistics 
(the  mean  and  variance  of  both  process  and  measurement  noise).  With  this  approach  used 
for  parameter  identification,  the  lower  bound  on  the  variances  of  parameter  estimates 
and  the  models  for  the  measurement  and  process  noise  disturbances  may  be  determined. 

The  majority  of  experience  gained  over  the  past  20  years  in  parameter  identification 
in  the  ballistic  range  and  wind  tunnel  has  been  accumulated  using  some  form  of  the 
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differential  correction  method  (see  Murphy,  Ref.  8,  and  Nicolaides,  Ref.  9,  for  example). 
The  maximum  likelihood  method  has  seen  limited  use  in  the  analysis  of  free-flight  aircraft 
test  data.  Thus,  even  though  the  maximum  likelihood  method  appears  to  hold  promise 
for  improved  data  analysis,  the  differential  correction  methods  as  used  in  the  wind  tunnel 
and  ballistic  range  are  stUl  of  primary  interest  due  to  their  relative  ease  of  application 
and  proven  reliability. 

The  differential  correction  method  as  used  by  Murphy  (Ref.  8),  Nicolaides  (Ref.  9), 
and  others,  is  essentially  a linear  or  slightly  nonlinear  analysis.  Under  suitable  linearizing 
assumptions,  the  angular  motion  equations  of  a basically  symmetric  model  may  be  solved 
in  closed  form.  The  coefficients  in  the  analytic  solution  which  are  functions  of  the  vehicle 
aerodynamics  are  determined  by  fitting  the  solution  to  the  data  using  the  method  of 
least  squares.  This  approach  can  be  used  only  to  extract  relatively  simple  aerodynamic 
coefficients  from  the  experimental  data.  Furthermore,  complex  motion  must  usually  be 
analyzed  in  small  time  segments,  and  thus  any  cumulative  nonlinearities  that  may  be  present 
in  the  data  are  destroyed. 

In  the  late  1960's  a different  approach  to  the  differential  corrections  method  was 
developed  that  did  not  rely  on  the  existence  of  analytic  solutions.  This  technique,  hereafter 
referred  to  as  the  Chapman-Kirk  method  (Ref.  5),  is  also  based  on  the  minimization  of 
the  least  squares  function.  The  partial  derivatives  as  required  in  the  method  of  least  squares 
are  determined  by  numerically  integrating  parametric  differential  equations  which  are 
derived  from  the  governing  equations  of  motion.  This  approach  can  satisfactorily  handle 
highly  nonlinear  aerodynamic  moments  and  forces,  but  it  requires  substantially  more 
computational  effort  since,  in  addition  to  the  integration  of  the  governing  equations,  a 
set  of  parametric  differential  equations  must  be  solved. 

The  discussion  of  output  error  methods  thus  far  has  dealt  with  parameter  identification 
techniques  that  are  deterministic  in  nature  in  that  the  aerodynamic  modeling  does  not 
account  for  measurement  or  process  noise.  All  ballistic  and  wind  tunnel  angular  motion 
data  contain  a finite  amount  of  measurement  and  process  noise,  however.  The  effect  of 
such  noise  usually  appears  in  the  higher  order  nonlinear  coefficients,  as  pointed  out  in 
Ref.  10,  in  the  form  of  abnormally  large  variance  and  bias.  The  extended  Kalman  filter 
method  (as  discussed  in  Ref.  1 1 ) provides  an  approach  in  which  measurement  noise  can 
• be  satisfactorily  modeled.  This  method  provides  direct  estimates  of  the  states  of  noisy 
systems  as  well  as  estimates  of  the  state  variable  uncertainties.  Even  though  the  extended 
Kalman  filter  method  provides  more  detailed  information  on  the  measurement  noise  of 
experimental  data,  it  provides  no  estimate  of  process  noise,  which  is  a serious  noise  source 
in  nonlinear  systems.  In  this  regard,  the  extended  Kalman  filter  method  has  no  distinct 
advantages  over  the  Chapman-Kirk  method  in  processing  real  3 DOF  data. 
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Thus,  the  Chapman-Kirk  method  is  used  to  develop  an  improved  3DOF  data  analysis 
capability  for  use  in  continuous  wind  tunnels.  Even  though  this  approach  does  not  yield 
detailed  information  on  measurement  noise,  as  does  the  extended  Kalman  filter  method, 
it  does  afford  sufficient  insight  into  the  effects  of  the  measurement  noise  on  the  extracted 
coefficients.  In  the  absence  of  process  (modeling)  noise,  it  may  be  shown  that  the  probable 
error  of  the  extracted  coefficients  for  random  measurement  noise,  using  the  Chapman-Kirk 
method,  is  equal  to  the  root-mean-square  (RMS)  deviation  from  the  true  value. 

3.0  MATHEMATICAL  DEVELOPMENT 

The  preceding  section  has  dealt  extensively  with  the  different  types  of  parameter 
identification  techniques  that  may  be  used  for  3DOF  data  analysis.  It  was  concluded  that 
the  Chapman-Kirk  method  as  described  in  Ref.  5 provided  the  best  compromise  in 
mathematical  complexity,  reliability,  and  versatility  of  the  methods  considered.  This  section 
presents  in  a concise  manner  the  mathematical  details  of  the  Chapman-Kirk  method  as 
applied  to  3 DOF  data  for  axisymmetric  bodies  of  revolution  taken  in  continuous  wind 
tunnels.  A similar  development  of  the  important  equations  for  models  with  mirror  plane 
aerodynamic  and  geometric  symmetry  is  given  in  Appendix  A. 

3.1  GOVERNING  EQUATIONS  FOR  3DOF  MOTION 

The  3 DOF  governing  equations  of  motion  are  derived  from  Newton's  second  law 
of  motion,  which  may  be  written  as 

M = iiii  ni 

y 

where  M is  the  external  moment  acting  at  the  mass  center  of  gravity  (eg)  equal  to  Mxi 
+ Mjj  + M/k,  and  h is  the  moment  of  momentum  referred  to  the  axis  system  as  shown 
in  Fig.  1.  For  a rigid  body,  the  moment  of  momentum  may  be  shown  (Ref.  12)  to  reduce 
to  Eq.  (2). 


h = + hyi  + h^k  = (pl^^  - “ '■^xz^  ' 

+ (-pixy  + ‘llyy  “ '‘^z^  j 
+ (_pI  - ql  + iT  ) k 

Here  p,  q,  and  r are  the  components  of  the  angular  velocity  vector  along  the  x,  y,  z 
axis  system  (see  Fig.  1),  respectively,  referenced  to  a space-fixed  axis  system,  and  l^x, 
lyy.  Iz/.>  Ixyi  Ixzi  lyz  moments  and  products  of  inertia  as  defined  in  the 

standard  way  by  Kolk  (Ref.  13)  and  others. 
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Ltilizing  Eq.  (2)  in  Eq.  (1).  the  governing  Euler  moment  equations  for  a plane-fixed 
axis  system  may  be  written  in  the  form  of  Eq.  (3),  where  the  caret  denotes  the  aeroballistic 


axis  system. 


A 


q h.^  - r h^, 

(3) 

pK  - 


The  dot  as  used  here  indicates  differentiation  of  the  scalar  quantity  with  respect  to  time, 

A A A 

and  p,  q,  and  r are  the  angular  velocities  of  the  plane-fixed  axis  system  relative  to  the 
space-fixed  axis.  For  the  case  in  which  the  body  axis  is  a principal  axis,  the  products 
of  inertia,  lx/.,  and  lyz,  are  identically  zero.  If  it  is  further  assumed  that  the  body 
is  axially  symmetric  about  the  x axis,  then  ly  = = 1. 

In  order  to  derive  the  Euler  moment  equations  in  terms  of  angular  measurements, 
the  Eulerian  angles  \p,  6,  and  0 are  defined.  Figure  2 illustrates  the  Euler  angles  8, 
and  4>  for  the  aeroballistic  axis  system.  Vectors  in  the  space-fixed  coordinate  system  may 
be  transformed  to  the  aeroballistic  axis  system  by  successive  rotations  about  the  Z axis 
through  the  yaw  angle,  \l/,  and  about  the  y axis  through  the  pitch  angle,  6. 

The  Euler  angles  \p,  6,  and  0 are  related  to  the  angular  velocities  through  Eq.  (4), 

p = (/}  - li/  sin  0 

q = 0 (4) 

r = ^ cos  0 


with  the  plane-fixed  angular  velocities  given  as  follows: 


P 


—ij/  sin  6 


T = Hi  cos  0 


(5) 


Inserting  Eqs.  (4)  and  (5)  into  the  governing  moment  equations  and  simplifying,  the 
governing  differential  equations  of  angular  motion  may  be  shown  to  reduce  to  the  form 
of  Eq.  (6). 


c 


II 


AEDC-TR -78-10 


c6  = ^'1  X - 6 cos  0 -I-  M ^ tan  0 - {2  - H)  t!j  6 tan  0 sin  6 

■ n 0 0 tan  0 

(6) 

0 = M^,  - tJ  - H)  sin  0 cos  0 - H w cos  6 


i!j  = M^/cos  0 I (2  - H)  0 0 tan  0 h-  R o 0/cos  0 

^ A A A 

Here  the  definitions  My  = My /I,  M,.  = jM,/!.  M^  = M^/Iv^  and  R = K/l  have  been  utilized. 
The  above  equations  describe  the  angular  motion  of  an  axially  symmetric  body  about 

A A A 

its  eg  caused  by  aerodynamic  moments  My,  and  M,,  in  the  aeroballistic  axis  system. 


The  solution  of  Eq.  (6)  depends  on  the  form  chosen  for  the  moments  as  functions 
of  the  dependent  variables.  In  order  to  provide  the  most  general  form  of  the  moment 
coefficients  for  axially  symmentric  bodies,  the  ideas  of  Tobak,  Ref.  14,  have  been  utilized. 
The  moments  are  written  as  functionals  of  the  parameters  that  directly  influence  the 
moments.  Thus,  it  is  assumed  that  the  moment  coefficients.  Cfi,  C,,, , and  Cn  are  functions 
of  the  velocity  ratios  v/U^  and  w/L'„  and  the  angular  velocities  p,  q,  and  r.  as  shown  in 
F-q.  (7). 


'1 X = C.^ia,  p-.  p , , r) 

My  = q^.Sd  (a.  p,  q,  r)  (7) 

Mj,  = Cj,  (a  , j3;  p,  q,  r) 


The  variables  a and  |3  are  the  velocity  ratios  w/U^  and  v/U„,  respectively,  as  depicted 
in  Fig.  1.  The  dynamic  pressure,  reference  area  (base  area),  and  reference  length  (base 
diameter)  are  denoted  by  q__,,  S,  and  d,  respectively.  By  expanding  the  moment  coefficients 
about  a and  |3  and  dropping  nonlinear  coefficients  in  a,  p,  q,  and  r,  the  coefficients 
of  Eq.  (7)  may  be  reduced  to 

(a  , /3  ; p , q , r)  = (a  . /3)  -i-  -A  ^ p C (a  , 13) 

h q (a,  /3)  -H  rC^^(a,  j8)  -i-  aC.^‘^(a,  f3)  + ^C^^(a. 

= C„,(a,^)  ^ [pC,^p(a,/3) 

(8) 

Cn(a,  p,  q,  rl  = C^^{a,^)  ^ 

+ qEn(,(a,^)  -I-  rC^^{a,f3)  + + /3C^^(a,/3)] 
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As  can  be  seen,  the  new  coefficients  as  defined  in  Eq.  (8)  are  functions  of  the  parameter's 
a and  j3  only.  Before  the  expansion  of  the  new  coefficients  about  a and  P equal  to  zero 
is  accomplished,  the  relations  describing  a and  P and  their  derivatives  as  functions  of  the 
Euler  angles  will  be  presented.  Using  this  information  will  result  in  simplifications  which 
ultimately  will  reduce  the  number  of  unknown  parameters.  Noting  that  the  only  component 
of  velocity  that  exists  in  a continuous  wind  tunnel  situation  in  the  space-fixed  coordinate 
system  is  -U„  along  the  X axis,  transformation  to  the  aeroballistic  axis  system  results 
in 

a = -Ji.  = cos  lb  sin  6 

I'. 


iS  = ^ 
U 


= -sin  ip 


(9) 


with  rates  of  change  of  a and  jS  given  by 

• • • 

a = -ip  sin  ^ sin  0 + 0 cos  ib  cos  0 

(10) 

j8  = -ip  cos  lif 

A comparison  of  Eqs.  (10)  and  (4)  makes  it  evident  that  for  small  pitch  and  yaw  angles, 
a is  approximately  equal  to  6.  This  of  course  assumes  that  6 and  fp  are  of  the  same 
order,  a valid  assumption  for  bodies  with  axial  symmetry  or  with  relatively  small 
asymmetries.  Thus,  the  rates  a and  P are  approximate  linear  combinations  of  the  angular 
rates  q and  r and  approach  exact  linear  combinations  in  the  limit  as  p and  6 approach 
zero.  Moreover,  in  Section  5.1  it  is  shown  that  for  reasonable  levels  of  measurement  noise, 
the  effects  of  the  rate  terms  a and  p cannot  be  satisfactorily  separated  from  the  angular 
velocities  q and  r.  Thus,  to  avoid  over  specification  of  the  aerodynamics,  it  is  assumed 
that  a = ej  r and  j3  = eir  where  ej  and  e2  are  constants.  This  reduces  the  coefficient 
expansion  as  given  in  Eq.  (8)  to  the  form 

Cj(a,/3;  p,  q,  r)  = Cj(a,/3)  + ^ j^pCjp(a,/3)  - qC;,.(a,/3)  i-  rCj^{a,/3)]  (11) 

where  i equals  (!,  m,  or  n. 

In  order  to  allow  for  a high  degree  of  versatility  in  the  coefficient  expansions  in 
terms  of  a and  P,  while  at  the  same  time  retaining  functional  simplicity,  Taylor  series 
expansions  are  made  about  a and  P equal  to  zero  for  each  of  the  coefficients  Q.  Qp, 
C,q,  and  Qr.  Terms  involving  the  magnitude  of  a and  |3  have  also  been  included  in  the 
expanded  coefficients.  Retaining  terms  cubic  and  lower  in  a and  p,  for  Cj(a,/3),  and 
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quadratic  and  lower  for  Qp,  Qq.  and  Qr,  the  expanded  form  of  the  coefficients  is  given 
in  Eq.  (12). 


C-ia, 

/3)  = 

il  • A -F  ^ ‘ 1 

lO  1 1 

a ■ 

r a ^ 

1 

il  1“ 

1 2 
. + ^i2  “ 

4 a'j2a  1 

a' 

“i3 

a 

i3 

a2 

ki 

+ ^14/3 

+ 

•*14  1 

^1 

+ 

a.. 

ID 

/32 

+ ^'isl 

1/81 

J8  4 

■ "iS 

/8® 

4 a 

16 

1/31/32  ^ 

a. 

1 

: 

+ 

a ' 

1 

7 1“ 

1/8 

4 a 

;7«i 

' 

i cl 

1 

iV  1 

a /3  1 4 a 

i8  “ 

2^ 

* 18 

«2l 

/81 

4- 

-.8  1 

a/3| 

a 1- 

“i9 

aP 

2 (■  e'jg 

1«1 

P-  4 

a 

i9 

fi  |a 

:/3| 

-/3)  = 

bjo  4 bj 

1“ 

•-  b 

il  1“ 

1 + 

*’i2' 

a2 

4 b'jjO 

l“l 

+ 

^3 

/3  4 

^ i3 

1/3 

1 + 

- b 

+ 

‘^3 

afi 

. b-slal/S  4 4 bV:  la^l 


(12) 


C„(a,i3)  = c 


iO 


[a  I 4 4 cj2aln| 


- 4 cjgl/SI  4 

. .\,\a\(i  * a.^l  cr^  'a/3' 


Cip(«./3) 


- ^iO  ~ 


f;,  a 4 f'l  |al  4 - fj,  a|a| 


- f.3^  - ■'  ■' 

- fV5«l/3-  4 r-\a^\ 


Here,  the  coefficients  a,,  b,.  Cj,  and  f,  are  independent  constants  that  are  functions  of 
the  free-stream  conditions  (such  as  Mach  number  and  Reynolds  number)  and  body 
geometry. 

Further  simplification  of  Eq.  (12)  is  not  possible  unless  flow-field  symmetry  arguments 
are  employed.  The  primary  justification  for  implementation  of  the  symmetry  conditions 
is  based  on  small  angular  motion  about  and  d equal  to  zero.  Since  the  maximum  angular 
displacement  of  the  3DOF  gas  bearing  system  is  10  deg,  flow-field  symmetry  that  is 
compatible  with  geometric  symmetry  is  implied.  Thus,  symmetry  arguments  as  given  by 
Maple  and  Synge,  Ref.  1 5,  are  used  to  reduce  the  number  of  unknown  constants  in  Eq. 
(12).  To  reduce  the  coefficients  as  given  in  Eq.  (12)  to  those  for  a body  symmetric  about 
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the  X axis  (see  Fig.  3),  two  covering  operations  are  required.  A covering  operation  for 
mirror  symmetry  about  the  x,  z plane  is  performed  first  and  is  followed  by  a covering 
operation  for  90-deg  rotational  symmetry  about  the  x axis.  The  mirror  symmetry  condition 
results  in 


-p.  q,-r)  = C^(a,^;p,  q , r) 

C„(a,-/3:-p,  q,-r)  = - C„  (a  , /3  ; p , q,  r)  (13) 

(a  , -jS  ; - p,  q,  -r)  ^ -C^  (a  , /3  ; p,  q,  r) 


Under  the  restrictions  of  Eq.  (13),  the  moment  coefficients  as  given  by  Eq.  (12)  reduce 
as  follows: 

C^{a,/3)  = - a^^,a  - a;,  a|  a ^ - a^jala] 


- + ®.n3  + “n, 41/31  - 


C.„q(a, 


,i31iQ2  H-  a 

';7«i/3'  - 

1/81 

+ 

**  mfi  \ CL  P ' 0. 

1/8^ 

+ ^.2«^  + 

b;„2a'a' 

+ Knsl^. 

- >^.4^^  - 

+ '^nVsl 

a^\ 

c,„>>/S)  = c;5ia|/3 

= ^..,3/3  + 

Ci(a,^l  = “io  + a, 4^  + + “ir,/3^  + ®.7“/3 


(14) 


-I  ajylal/S  + a.ga^^  + a'jg/Sla^l 
C.^(a,^)  = bjg^  -I-  bj^/SI/Sl  1 bjga/S  + b^glaj/S 

Cir(a./3)  = CjQ  ^ c ija  + c i,  >1  + c.^a^  + c\^a  ,a\ 

+ + ^14/3'  - c'isal^'  - 

Cip(a,j8)  fjo  V ri,a  - filial  1 ■-  fjaalal 

- f;3l/3  •-  f'i'5«liS|  + fi5l“/31 
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Here,  i = n,  C.  The  coefficient  expansions  as  given  in  Eq.  (14)  are  valid  for  asymmetric 
configurations  with  the  symmetry  plane  of  the  body  in  the  x,  z plane  of  the  body-fixed 
axis  system.  A 90-deg  rotational  symmetry  requirement  will  result  in  a four-fold  symmetry 
of  the  body.  The  four-fold  symmetry  coupled  with  mirror  plane  symmetry  is  sufficient 
for  the  specification  of  axial  symmetry.  A covering  operation  performed  on  a 90-deg 
rotation  of  the  axis  system  results  in  the  coefficient  constraints  as  given  in  Eq.  (15). 


C„(^,  -a)  = C^(a, 

()3  , -a)  = C|  (a  , /3) 
-C/qC/3  , -a)  = /3) 

C^p(/3,  -a)  = C^p(a,/3) 

= C,p(a,^) 


Employing  these  contraints  with  the  additional  restriction  that,  for  zero  spin  rate,  an  axially 
symmetric  body  will  have  no  moments  in  the  plane  normal  to  the  total  angle-of-attack 
plane,  the  aerodynamic  coefficient  expansions  as  shown  in  Eq,  (16)  will  result. 

C^Q(a,0)/a  = aj  + + [ /3  j ) + ^3  (a  ^ + )S^)  + a''\aP\ 

- €5/32  = bo  + bjdal  + 1)31)  + b2(a2  + + b'''|a)8| 

= ‘^3  + ‘‘sl®'!  + **41^1 

-Cno(a. /il)//3  = 

C„q(«./3)  = 

C„p(a,)8)/a  = dg  + d:i)3|  + d^'a] 

Cjp(a,)3)  = fo  + f'ldal  + li8|)  + f2(a2  ^ 

C^r(a,)3)/a  = g3  - g'sl^li  + g4la' 

Ci(.(a,/3)//3  = g3  + gglal  + g4|/3, 
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The  above  expressions  as  used  in  Eq.  (11)  yield  the  required  moment  coefficients. 

It  is  important  to  note  here  that  the  expansions  as  given  in  Eq.  (16)  are  not  unique 
in  describing  axially  symmetric  body  aerodynamics.  Indeed,  other  expansions,  such  as  those 
in  terms  of  the  total  angle  of  attack  as  used  by  Whyte  (Ref.  16)  and  others,  are  also 
valid.  The  "correct"  expansion  of  coefficients  describing  the  aerodynamics  of  a body, 
in  the  parameter  identification  sense,  is  that  which  yields  the  "best"  fit  to  the 
experimental  data  and  does  not  violate  established  symmetry  constraints.  With  the 
establishment  of  the  required  aerodynamic  coefficients,  the  governing  equations  of  motion, 
Eq.  (6),  may  be  solved. 

3.2  FORMULATION  OF  PARAMETER  EXTRACTION  METHOD 

In  order  to  determine  the  coefficients  describing  the  aerodynamic  model  as  given 
in  Eq.  (16)  for  axial  symmetry  or  in  Eq.  (14)  for  mirror  plane  symmetry,  a parameter 
extraction  technique  must  be  used.  The  3 DOF  gas  bearing  as  described  in  Ref.  1 yields 
angular  position  measurements  in  yaw,  pitch,  and  roll  which  provide  the  experimental 
data  necessary  to  extract  the  unknown  parameters.  In  Section  1.0  the  Chapman-Kirk 
method  was  selected  as  the  parameter  identification  technique  to  be  employed  here.  This 
technique  is  an  output  error  method  and,  as  such,  uses  the  method  of  least  squares  to 
minimize  the  cumulative  error  between  observed  and  calculated  motion. 

Normally  when  one  is  fitting  data  with  a least  squares  method,  the  function  to  be 
fitted  is  given  with  constant  coefficients  that  appear  linearly  in  the  function.  This  results 
in  a set  of  simultaneous  linear  equations  that  may  be  solved  in  a noniterative  fashion 
for  the  unknown  coefficients.  For  nonlinear  functions,  the  coefficients  must  be  evaluated 
by  the  differential  corrections  method.  This  technique  provides  corrections  to  improve 
the  estimates  of  the  unknown  parameters  and,  as  a result,  is  an  iterative  process. 

Since  the  3 DOF  governing  equations  as  given  in  Eq.  (6)  are  nonlinear  and  the  unknown 
parameters  as  given  in  Eq.  (16)  are  complicated  functions  of  the  Euler  angles,  4>,  0,  and 
\}j.  the  differential  corrections  method  will  be  used.  Consider  the  problem  of  fitting  the 
3DOF  motion  of  a model  to  known  functions  of  the  unknown  parameters.  If  the  variables 
0m  £■  and  are  defined  as  the  experimental  values  of  the  roll,  pitch,  and  yaw 

positions  at  each  time  point  residuals  (ug,  vg,  wg)  yielding  the  difference  between  the 
measured  and  calculated  data  may  be  expressed  by  Eq.  (17). 

(17) 

= 0/  - 
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Here,  pfi,  08,  and  i//g  are  the  calculated  values  of  the  roll,  pitch,  and  yaw,  respectively, 
which  are  functions  of  a set  of  independent  constants,  a,.  These  constants  are  the 
aerodynamic  coefficients  to  be  determined  as  shown  in  bq.  (16).  In  order  to  determine 
the  change  in  the  coefficients,  required  to  reduce  the  residuals  in  the  least  squares 
sense,  a Taylor  series  expansion  of  the  calculated  solution  is  made.  Expanding  Os, 
and  i//S,  and  discarding  second  and  higher  order  terms  and  using  the  definitions  as  given 
in  Eq.  (17),  one  obtains  Eq,  (18). 


N 

r 


I 


N 

c 


(18) 


+ W 


mi 


^io 


+ T Aa.  ^ 
^ ' Oa 

1=1  * 


The  functions  ^Co.  and  i/^So  correspond  to  the  values  of  0S,  and  i//g,  respectively, 

evaluated  at  each  time  point  for  the  given  constants  aj.  and  N'c  is  the  total  number  of 
parameters  a,  to  be  determined.  The  best  fit  of  ipSi  0g,  <ind  i//S  to  the  experimental  data 
is  established  by  minimizing  the  sum  of  the  squares  of  the  residuals  for  all  time  points 
considered.  Thus, 


is  subject  to  the  minimization  condition  as  given  in  Eq.  (20). 

N 


0A.. 


2 ^“1  + '■j  + = 0 


(20) 


i = i 


where 


^ ' P2ii-  Psu- 
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N is  the  total  number  of  time  points  considered,  and 

^1/  “ ~ ■ ^2i  = ~ ■ J^3i  ‘ '^-ni  “ io 

Performing  the  indicated  partial  differentiation  in  Kq.  (20)  and  combining  terms  results 
in 


C,.  All.  ^ I) 


(21) 


where  C,j  is  a symmetric  (N^  x Nj.)  matrix  of  parametric  influence  coefficients  given  by 

N 

<'-i  - EC’iuPiji  - 1 - PajPsi.)  ™ 

/ = 1 

and  Dj  is  a residual  influence  matrix  with  the  form 

N 

«j  - ZO'^Pix-  I’aiPsji)  (23) 

i-1 

Since  Eq.  (21)  is  linear  in  Aa,,  it  may  be  inverted,  yielding  the  required  correction  to 
the  coefficients,  a,.  Thus,  given  a set  of  initial  estimates  for  the  independent  constants, 
determination  of  Aa,  will  result  in  corrected  values  of  a,  which  will  cause  the  computed 
motion  to  more  closely  approach  the  experimental  data.  In  the  limit  of  no  process  or 
measurement  noise,  iteration  on  a,  will  provide  the  best  fit  in  the  least  squares  sense 
of  the  computed  to  the  experimental  motion. 

Implied  in  this  development  are  two  important  facts  that  can  seriously  disrupt  or 
destroy  the  convergence  of  a;  to  the  "best"  values.  If  measurement  noise  exists,  the 
definition  of  the  residuals,  Eq.  (19),  assumes  that  the  dependent  variables  0,  6,  and 
have  the  same  relative  uncertainty.  For  practical  applications,  the  relative  uncertainty  of 
the  measurements  <i>.  6,  and  i//  may  differ  by  as  much  as  an  order  of  magnitude.  This 
effect  can  be  accounted  for  in  the  parametric  influence  coefficients  by  redefining  Cy  as 


and  Dj  as 


p p 
‘ 2ii  ' 2ii 

^0 


4 ) 


^6 


4 ) 


(24) 


(25) 
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where  a^.  0$,  and  are  the  standard  deviations  of  the  relative  uncertainty  in  the  roll, 
pitch,  and  yaw  measurements,  respectively  (see  Ref.  17).  Also  implied  in  Eq.  (18)  is  the 
assumption  that  the  coefficients  (a,)  used  to  generate  0£o,  0£o,  and  Mo  are  close  to 
the  best  values.  If  a significant  difference  exists  between  the  initial  value  of  the  coefficients 
and  the  actual  or  best  values,  the  expansions  as  given  in  Eq.  (18)  are  not  valid  and  will 
result  in  divergence  of  the  corrections,  Aa,.  or  convergence  on  invalid  solutions.  The 
closeness  of  the  initial  values  of  the  coefficients  to  the  best  values  for  correct  convergence 
is  a function  of  the  character  of  the  motion  fitted. 


As  stated  previously,  in  the  absence  of  measurement  or  process  noise,  correct 
application  of  the  differential  corrections  method  will  yield  parameters  with  no  relative 
error.  In  the  presence  of  these  noise  sources,  however,  bias  may  be  introduced  into  the 
converged  parameters.  According  to  Ref.  1 8,  for  statistically  random  measurement  noise 
in  the  data,  it  may  be  shown  that  the  root-mean-square  (RMS)  deviation  betw-een  the 
experimental  data  and  converged  solution  is  given  by 


E = 


1)  - N n 

r 


with  the  corresponding  RMS  error  in  each  of  the  coefficients,  Saj,  equal  to 


(26) 


where  Bi  represents  the  diagonal  components  of  the  inverse  matrix  of  Qj.  Equations  (26) 
and  (27)  provide  mathematically  rigorous  estimates  of  the  measurement-noise-generated 
R.MS  deviation  in  the  coefficients  and  solution.  They  do  not.  however,  give  reliable 
estimates  of  the  RMS  deviation  or  error  for  measurement  noise  that  is  not  statistically 
random,  noise  which  may  result  when ' data  time  samples  are  too  small,  or  for  process 
noise  caused  by  inappropriate  aerodynamic  modeling. 

The  differential  corrections  method  just  presented  requires  the  value  of  each  of  the 
partial  derivatives,  as  illustrated  in  Eq,  (18),  at  each  time  point.  Previous  users  of  the 
differential  corrections  method  such  as  Eikenberry  (Ref.  2)  relied  on  the  use  of  analytic 
solutions  for  the  angular  motion  to  generate  the  required  partial  derivatives.  Since  the 
equations  of  motion  as  used  here  are  highly  nonlinear  and  cannot  be  solved  in  closed 
form,  the  values  of  the  influence  coefficients  (di^/daj,  90/9ai,  and  9 \!//9ai)  are  determined 
by  numerically  integrating  a set  of  parametric  differential  equations.  These  equations  are 
derived  by  differentiating  the  governing  equations  of  motion  with  respect  to  each 
parameter,  a,,  and  inverting  the  order  of  differentiation.  Performing  the  indicated  operations 
on  Eq.  (6)  results  in  the  following: 
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Pjj  = M|j  + 0 Pgj  COS  0 + i!/  Pjj  cos  0 - 0 P2i  sin  9 

+ Mgi  lan  6 + P2i  sec.^  0 + (2  - K)  (0P3,  -i-  (JfPjj)  tnn  6 sin  0 

+ (2  - R)  ^ffPgj  sin  0 (1  - sec“  0)  : R 0 P|  j lan  0 

+ R0P2jtanO  I-  R(i0P,jSec^0 

P,.  = Mgi  - 2 (1  - R)  0P3i  sin  0 cos  0 - (1  - R)  ijj~  Pg;  (1-2  sin^  0)  (28) 

- R P^  j cos  0 - R (iPgj  cos  0 h R ^ 0 P2,  sin  0 

Pji  = iVIgj  sec  0 + My,  1^2;  tan  0 sec  0 + (2  - R)  tan  0 (i,Z'P2,  *•  0 ^31) 

(2  - R)  ^0P2;  sec^  0 -I-  R 0 P[j  see  0 + R ©0  tan  0 sec  0 P2[ 

+ R <i  P,;  sec  0 

aM^  aMy  awz 

where  M] ; = , M2i  = and  M3i  = are  evaluated  in  a similar  manner  through 

dai  aaj  aaj 

the  moment  expansions.  The  above  equations  are  second  order  and,  thus,  require  two 
boundary  conditions  for  the  evaluation  of  each  influence  coefficient.  These  are  obtained 
from  the  boundary  conditions  on  the  angular  motion  given  by 


d(0)  = ttjj  _3 
(• 

(0)  = Qjij  _ ^ 

c 

0 (0)  = a _ 3 

0(0)  = aj^.  _2 

c 

1^(0)  = 


(29) 


^(0)  = ajj__ 

If  one  differentiates  the  boundary  conditions  given  in  Eq.  (29),  the  boundary  conditions 
on  the  governing  equations  for  the  influence  coefficients  may  be  derived. 
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Pl.^0)  = Sn  _5,i 

r 

P,  ;(0)  = Sjj  _ i 
c 

PjjCO)  = _ 3 i 

c 

PgiCO)  = 5p^.  _2_i 

C 

Psi^o^  = 5n  -1., 

c 

P3i(0)  = Sj^,  j 

c 


(30) 


Here  Sjj  is  the  Kroneker-Delta  function  with  8,j  =1  for  i = j and  S,j  = 0 for  i j. 
In  order  to  utilize  the  differential  corrections  method  to  extract  aerodynamic  coefficients 
of  interest,  one  must  solve  the  governing  angular  motion  equations  and  parametric 
differential  equations  simultaneously.  With  the  angular  motion  and  influence  coefficient 
solutions  determined  for  a given  set  of  initial  coefficient  values,  corrections  to  the 
coefficients  are  evaluated  through  Eq.  (21).  For  example,  if  there  are  thirty  coefficients 
to  be  identified  [which  include,  generally,  both  the  boundary  conditions,  Eq.  (20),  and 
the  aerodynamic  coefficients,  Eq.  (I6)J,  a maximum  of  186  differential  equations  must 
be  solved  simultaneously. 


4.0  DESCRIPTION  OF  COMPUTER  PROGRAM 

The  computer  program  as  described  here  extracts  aerodynamic  coefficients  from  3 DOF 
data  defined  by  the  measured  yaw.  pitch,  and  roll  angles  as  a function  of  time.  The  code 
is  based  on  the  Chapman-Kirk  technique  as  detailed  in  Section  3.0.  Both  symmetric  and 
asymmetric  versions  of  the  program  have  been  developed  and  are  listed  in  Tables  1 and 

2. 

The  symmetric  version  is  written  in  double  precision  Fortran  IV  for  use  on  the  IBM 
370/165  computer.  In  the  form  listed  here,  core  requirements  are  164  K bytes  using  the 
G compiler.  As  many  as  30  coefficients  may  be  extracted  from  3DOF  data  with  actual 
computation  time  proportioned  to  t as  shown  in  Eq.  (31). 

r « n(.\'„  + 1)  (31) 

Here  n is  the  number  of  iterations,  Ne  is  the  total  number  of  constants  varied,  tmax 
is  the  total  integration  time,  and  At  is  the  integration  time  interval. 

The  program  is  written  in  a logical  sequence  so  that  changes  in  internal  program 
structure  such  as  alteration  of  the  functional  form  of  the  moment  coefficient  expansions 
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can  be  made  easily.  A flow  chart  showing  the  important  logic  and  branching  points  is 
provided  in  Fig.  4.  Table  3 contains  a brief  description  of  each  of  the  subroutines  employed 
in  the  code. 

The  asymmetric  version  of  the  code  is  functionally  the  same  as  the  symmetric  version. 
Only  the  core  requirements,  coefficient  arrays  in  the  main  program,  and  subroutines  directly 
involving  the  kinematic  or  dynamic  equations  for  the  3DOF  motion  are  altered.  Core 
requirement  for  the  asymmetric  program  is  266  K bytes  when  the  G compiler  is  used. 
In  addition,  up  to  72  coefficients  can  be  extracted  from  3 DOF  data  using  this  version. 

In  both  programs,  integration  of  the  kinematic,  dynamic,  and  parametric  differential 
equations  is  performed  with  a fourth-order  Runge-Kutta,  fixed  time  step  scheme.  The 
coefficients  to  be  extracted  from  the  data  are  reiterated  as  outlined  in  Section  3.2  until 
convergence  is  achieved.  Convergence  is  assumed  when  consecutive  iterations  yield  a change 
in  the  probable  error  of  the  fit  less  than  a specified  error  bound. 

4.1  INPUT 

The  input  information  necessary  to  execute  either  version  of  the  coefficient  extraction 
program  is  very  similar.  The  actual  3 DOF  experimental  data  to  be  analyzed  is  input  to 
both  programs  through  a direct  access  input  file.  Data  are  stored  in  this  file  in  1,003 
word  records,  in  which  the  first  three  words  define  the  data  group  number,  data  block 
number,  and  data  type  indicator,  in  that  order,  and  are  integers.  The  data  type  indicator 
defines  the  variable  being  accessed  and  may  take  on  values  from  one  to  four  corresponding 
to  arrays  of  time,  0 versus  time,  d versus  time,  and  ^ versus  time,  respectively,  all  of 
which  are  single  precision  variables. 

The  initial  values  of  the  coefficients  and  the  parameters  providing  branching 
instructions  and  body  and  flow-field  information  are  input  using  namelist  statements.  No 
default  value  is  provided  in  the  program  for  required  variables  not  defined  in  the  initial 
namelist  statement.  The  namelist  input  format  was  used  to  provide  versatility  in 
constructing  sequential  runs  by  updating  prior  information.  Table  4 gives  the  namelist 
variables,  their  meaning,  and  their  nominal  values  needed  for  execution  of  either  the 
symmetric  or  asymmetric  coefficient  extraction  program. 

4.2  OUTPUT 

Output  variables  as  used  in  both  the  symmetric  and  asymmetric  programs  are  similar 
and  are  compatible  with  input  variable  definitions.  Only  portions  of  the  output  data  from 
a typical  run  which  are  not  self  explanatory  will  be  discussed. 
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During  normal  data  reduction  using  either  program,  each  iteration  of  the  parameters 
(aerodynamic  coefficients  and  boundary  conditions)  is  listed  with  the  corresponding  change 
in  the  parameter  (DELTA  CON)  according  to  the  theory  of  Section  3.2.  The  variables 
defined  as  the  new  and  old  probable  error  are  the  previous  and  present  values  of  the 
RMS  deviation  between  the  experimental  data  and  converged  solution,  respectively,  as 
defined  in  Eq.  (26).  The  sum  of  the  residuals  is  defined""as  follows: 


Sum  of  Residuals  ''  (32) 

The  iterations  are  continued  until  convergence  of  the  solution,  as  discussed  in  Section 
4.0,  is  achieved  or  until  the  input  value  of  the  maximum  number  of  iterations  is  reached. 
The  final  values  of  the  parameters  are  then  listed  along  with  corresponding  values  of  the 
RMS  error  in  each  coefficient  as  defined  in  Eqs.  (26)  and  (27).  A comparison  of  the 
fit  of  the  theoretical  solution  (denoted  by  FIT)  to  the  experimental  data  (denoted  by 
EXP)  is  also  shown.  The  theoretically  generated  solution  of  ^ and  0 (denoted  by  PSI  and 
THETA,  respectively)  and  the  difference  between  the  theoretical  and  experimental  values  of 
\l)  and  0 i\jj  - \!/m  = PSI  - PSIM,  0 - 0m  = THETA  - THETAM)  are  stored  for  later  use  as  data 
for  graphical  output. 


5.0  RESULTS  OF  VALIDATION  TESTS 

Verification  of  the  improved  parameter  identification  technique  as  presented  in 
Sections  3.0  and  4.0  requires,  in  addition  to  checks  on  the  program  using 
computer-generated  motion,  an  assessment  of  the  measurement  and  process  noise  inherent 
in  the  gas  bearing  and  in  typical  tunnel  test  data.  This  is  necessary  since  the  "best" 
coefficients,  as  produced  by  output  error  methods  such  as  the  one  discussed  here,  are 
directly  influenced  by  the  cumulative  effects  of  modeling  errors  and  measurement  noise. 
Thus,  this  section  addresses  three  areas  concerned  with  validation  of  the  improved 
coefficient  extraction  technique:  numerically-generated  3 DOF  motion,  including  Gaussian 
noise,  is  employed  to  simulate  measurement  noise  effects  on  the  probable  error  in  extracted 
coefficients;  a bench  test  is  used  to  assess  the  accuracy  of  the  3DOF  measurement  system; 
and  3DOF  data  taken  in  Tunnel  A are  used  to  indicate  some  of  the  capabilities  of  the 
parameter  identification  technique. 

5.1  COMPUTER-GENERATED  MOTION 

In  order  to  arrive  at  the  final  form  of  the  aerodynamic  coefficients  as  given  in  Eq. 
(14)  for  the  asymmetric  and  Eq.  (16)  for  the  symmetric  analysis,  the  a and  |3  coefficients 
were  linearly  combined  with  the  q and  r rate  coefficients,  respectively.  This  was  a direct 
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result  of  similarity  of  the  a,  q and  jS,  r relations  [(see  Eqs.  (4)  and  (10)]  in  terms  of 
the  Euler  angles.  In  order  to  verify  this  hypothesis,  3DOF  motion  was  generated  botli 
with  and  without  noise  with  the  a and  j3  coefficients  included.  The  initial  conditions  and 
parameters  employed  to  specify  the  motion  are  listed  below. 

S = 0.390  ft2,  = 0.0372  slugs-ft2,  1 = 0.244  slugs-ft2 

Cma  = -0.236  1/rad,  Cmq  = -1.762  1/rad,  Cmi  = -0.881  1/rad 

0o  = 0,  00  ” 641  deg/sec,  do  = 5.18  deg,  0„  = -26.356  deg/sec 

\l/o  = -2.589  deg,  i0o  = -7.77  deg/sec 

= 4,774  ft/sec,  q„  = 150.4  lbf/ft2 

The  symmetric  version  of  the  parameter  identification  program  was  used  to  extract  the 
variables  listed  above  from  the  generated  data.  The  parameters  q and  Cm  i were  both 
extracted  with  nominal  values  of  -1.74  per  radian  and  -0.90  per  radian,  respectively,  and 
RMS  deviations  of  ±0.05  per  radian.  This  deviation  is  a measure  of  the  character  of  data 
fitted  and  the  numerical  accuracy  of  the  program.  Gaussian  noise  with  one  standard 
deviation  equal  to  0.01  deg  was  then  added  to  the  generated  values  of  <j>,  d,  and  0.  With 
noise  added,  Cmq  and  Cm  a were  both  extracted  from  the  data  and  yielded  values  of 
-5.87  per  radian  and  3.2  per  radian,  respectively,  with  R.MS  deviations  of  ±18  per  radian. 
The  noise  added  is  representative  of  the  levels  of  measurement  noise  inherent  in  the3DOF 
measurement  system.  Moreover,  the  motion  examined  contained  only  linear  aerodynamics. 
This  should  have  resulted  in  the  smallest  RMS  deviations  of  the  extracted  coefficients 
as  influenced  by  modeling.  Thus,  based  on  the  results  of  this  example,  it  is  not  feasible 
to  separate  the  effects  of  Cmq  and  Cmd  typical  3DOF  measurements  where  1-a 

measurement  errors  are  0.01  deg  or  higher. 

The  definitions  of  a and  jS  can  also  play  a role  in  the  extraction  of  the  "best" 
coefficients  for  a given  set  of  3DOF  measurements.  In  Tobak's  developments,  Ref.  14, 
a and  jS  were  defined  as 


a 


V 


V 


(33) 


Much  of  the  past  work,  however,  has  employed  a and  in  terms  of  the  angle  of  attack 
and  sideslip  leading  to . 


«a  = tan  > ^ , /3„ 


tan 


i] 


(34) 
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Written  in  terms  of  the  Euler  angles  in  the  aeroballistic  axis  system,  ay,  jSy  may  be  related 
to  Oa,  |3a,  resulting  in 


a = ct 

V 

1 

(■  i 

cos  ip 

6 

COS^  iff 

Pa  = 

\ 

1 

— + L 

1 cos  0 

COS 

xfj  3 

/3 


0 iff 


(35) 


Thus,  for  small  values  of  6 and  Oa  « Ov  and  Pa  ^ Pw-  Over  a large  number  of  cycles 
of  data,  the  cumulative  effect  of  these  small  differences  will  appear  as  higher  order  terms 
in  the  coefficient  expansions  [compare  Eqs.  (35)  and  (16),  for  example].  Both  definitions 
of  a and  P as  given  in  Eqs.  (33)  and  (34)  have  been  incorporated  in  the  parameter 
identification  programs  as  developed  here.  Comparisons  of  generated  data  using  the  two 
definitions  of  a and  P were  made  which  confirmed  the  results  derived  from  Eq.  (35). 

The  accuracy  and  validity  of  generated  3DOF  motion  from  the  programs  as  listed  in 
Tables  1 and  2 were  examined  througli  comparison  with  results  from  an  existing  6DOF 
code,  Ref.  19.  Calculations  were  made  using  the  6DOF  code  and  the  asymmetric  program 
as  described  in  this  report  for  two  sets  of  initial  body  and  flow-field  conditions  as  listed 
below. 

CHECK  CASE  1 
S = 0.3898  ft2,  d = 0.7038  ft 

\y  = h = 0.2444  slugs-ft2,  Ix  = 0.03723  slugs  -ft^,  l^z  = 0 
= 4,774  ft/sec,  q„  = 364.66  lbf/ft2 
Cm  a = -Cn(j  = -0.2355  1/rad,  C„q  = Cnr  = -2.6434  1/rad 
Cgp  = -0.00745  1/rad, 

do  = 5.18  deg,  \l/o  = -2.589  deg,  0o  = 0,  Po  = 641.7  deg/sec 
qo  = -26.356  deg/sec,  Tq  = -7.735  deg/sec 

CHECK  CASE  II 


S = 1.228  ft2,  d = 1.25  ft 

ly  = = 700  s!ugs-ft2.  = 20  slugs-ft2,  = -100  slugs-ft2 
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= 1,500  ft/scc,  = 675  Ibf/ft^ 

Cm  a = -Cns  = 16.7132  !/rad 

do  = 5 deg,  \l/o  = 5 deg,  = 0.  Po  = do  = fo  = 0 


Typical  results  of  the  calculations  are  shown  below  and  indicate  excellent  agreement 
between  the  two  codes  in  the  computed  values  of  0,  0,  and 


Results: 

Check  Case  I 

t, 

e. 

0, 

sec 

deg 

deg 

deg 

Present  3DOF  Code 

30 

-4.270 

0.702 

-23.798 

6DOF  Code  (Ref.  19) 

30 

-4.2704 

0.7001 

-23.762 

Results: 

Check  Case  11 

t, 

0. 

sec 

deg 

deg 

deg 

Present  3 DOF  Code 

19.90 

-0.429 

5.388 

25.766 

6DOF  Code  (Ref.  19) 

19.90 

-0.4289 

5.3883 

25.7659 

To  attain  the  best  possible  coefficients  from  3DOF  data  with  measurement  noise, 
the  quantity  and  character  of  the  data  sampled  are  of  fundamental  importance.  A less 
than  optimum  data  sample  will  yield  extracted  coefficients  with  high  RMS  deviations  or, 
in  some  cases,  lack  of  convergence  of  the  solutions.  The  extraction  of  optimum  coefficients 
was  studied  using  the  generated  motion  as  discussed  previously  with  0.05-deg,  1-a  Gaussian 
noise  added  to  the  data  and  = 0,  Cmq  = -2.64  per  radian  where  a is  one-standard 
deviation.  Table  5 shows  the  results  of  this  study  where  the  RMS  deviation  in  percent 
of  Cmq  is  shown  as  a function  of  the  number  of  data  points  per  cycle  of  data  and  the 
total  number  of  cycles  of  data  analyzed.  From  this  table  it  may  be  generally  concluded 
that  accuracy  in  the  extracted  value  of  Cmq  is  improved  faster  by  increasing  the  data 
sample  size  than  by  increasing  the  number  of  points  per  cycle.  This  same  trend  and  the 
rate  of  improvement  of  the  extracted  parameter  was  also  found  to  hold  true  for  Cm  a 
and,  moreover,  for  the  nonlinear  expansions  of  Cmti  and  Cmq  [see  Eq.  (16)]. 

The  effect  of  the  magnitude  of  the  measurement  noise  on  the  extracted  coefficients 
was  examined  by  generating  noisy  3 DOF  data  for  the  same  conditions  as  previously  used 
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but  with  Cm  as  2 = S.Sl/rad^,  Cmq  = -4.0  per  radian,  and  Cmqs^  = 98.5/rad3.  The 
symmetric  version  of  the  parameter  identification  program  was  used  and  yielded  the  results 
shown  in  Table  6.  Here  the  extracted  parameters  and  their  RMS  deviation  are  shown 
as  a function  of  the  measurement  noise  added  to  the  generated  values  of  0,  8,  and  0. 
The  generated  data  analyzed  contained  40  cycles  of  data  with  20  points  per  cycle.  Initial 
parameter  estimates  were  set  at  Cm  a “ -0.2  per  radian  and  Cmq  = Cmqs^  = Cmas^ 
= 0.  Typical  plots  of  calculated  6 and  0 values  versus  time  for  a measurement  noise 
level  of  a = 0.05  deg  are  shown  in  Fig.  5.  Residuals  showing  the  difference  between 
the  generated  data  with  noise  and  the  calculated  values  of  6 and  0 and  graphs  indicating 
the  fit  of  the  computed  solution  to  the  generated  data  are  also  provided. 

Two  important  conclusions  may  be  drawn  from  the  magnitude  and  variation  of  the 
RMS  deviation  of  each  of  the  coefficients  in  Table  6.  It  is  apparent  that  the  magnitude 
of  the  RMS  deviation  is  approximately  a linear  function  of  the  measurement  noise  level 
denoted  by  a.  In  addition,  the  . magnitude  of  the  deviation  appears  to  be  a good  estimate 
of  the  actual  error  in  the  values  of  the  extracted  coefficients.  This  result  is,  in  fact, 
verification  of  those  statements  made  in  Section  3.2  concerning  the  theoretical  meaning 
of  6ai  [Eq.  (27)]. 

5.2  BENCH  TEST  DATA 

The  previous  section  discussed  the  effects  of  measurement  noise  on  the  extracted 
coefficients.  The  noise  source  considered  was  Gaussian  and  random  in  nature. 

When  3DOF  motion  including  Gaussian  noise  is  viewed  in  the  frequency  plane,  the 
noise  contributes  no  distinct  frequencies  or  regions  of  high  spectral  density.  In  most  cases, 
however,  the  effect  of  modeling  noise  on  the  extracted  coefficients  is  more  pronounced 
than  the  effect  of  measurement  noise.  Most  modeling  sources  occur  over  finite  spectrums 
with  median  frequencies  occuring  near  aerodynamic  frequencies  of  interest.  Moreover,  in 
such  regions  modeling  noise  may  have  high  power  levels  which  may  be  the  same  order 
of  magnitude  as  noise  caused  purely  by  aerodynamics.  Thus,  a critical  assessment  of  the 
important  modeling  errors  present  in  typical  3DOF  data  must  be  made  to  gain  the  optimum 
capabilities  of  the  improved  parameter  identification  technique. 

Some  of  the  important  categories  of  modeling  noise  are  1)  incorrect  aerodynamic 
modeling  of  the  specific  configuration,  2)  modeling  errors  introduced  through  the  effects 
of  tunnel  flow  nonuniformities,  and  3)  errors  associated  with  nonlinearities  and  coupling 
effects  caused  by  the  gas  bearing.  This  section  will  consider  modeling  errors  due  to  the 
gas  bearing  only. 
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The  VKF  spherical  gas  bearing  as  described  in  Ref.  1 provides  a nearly  frictionless 
pivot  in  roll,  pitch  and  yaw,  and  thus  is  ideal  for  3DOF  motion  simulation.  Model 
oscillations  of  up  to  10  deg  in  the  combined  pitch,  yaw  plane  and  unlimited  roll  orientation 
are  possible  using  this  mechanism.  Figs.  6 and  7 are  photographs  of  the  gas  bearing 
assembly.  Details  of  the  balance  mechanism  and  the  angular  measurement  system  are  shown 
in  Fig.  8.  Pitch  and  yaw  angular  displacements  arc  measured  using  an  eccentric  as  shown 
in  Fig.  7 positioned  on  the  outer  race  of  the  bearing  and  two  pairs  of  orthogonal  E-corcs 
fixed  on  the  bearing  instrument  ring.  Motion  in  pitch  or  yaw  causes  a change  in  the 
distance  between  the  eccentric  and  the  E-corc  surfaces.  The  accompanying  change  of  the 
reluctance  of  the  E-cores  provides  an  analog  signal  proportional  to  angular  displacement. 
The  measurenaent  of  roll  position  is  made  using  light-emitting  diodes  and  alternating 
reflective  and  nonreflective  lines  affixed  to  the  eccentric  as  shown  in  Figs.  7 and  8. 

The  measurement  of  angular  position  in  pitch  and  yaw  using  the  gas  bearing  is  affected 
by  two  sources  of  error.  During  calibration,  adjustments  to  the  bearing  are  made  so  that 
the  E-core  output  varies  approximately  linearly  with  the  actual  angular  displacement. 
Calibration  factors  are  then  established  for  both  the  pitch  and  yaw  measurements  using 
a linear  least  squares  curve  fit.  Although  this  approach  provides  a very  simple  and  fast 
technique  for  reducing  raw  data  into  angular  displacements,  it  can  introduce  error  into 
the  resulting  measurements.  Nonlincarities  on  the  order  of  0.02  deg  or  less  for  pitch  and 
yaw'  measurements  of  up  to  ±5  deg  are  typical  for  bearing  calibration.  For  pitch  or  yaw 
measurements  approaching  8 deg,  the  nonlinearity  generally  increases  in  magnitude, 
reaching  values  as  high  as  0.1  deg. 

An  additional  source  of  error  is  the  interaction  of  roll  position  on  the  pitch  and 
yaw  measurements.  This  effect  is  most  probably  a result  of  the  machining  tolerances  of 
the  eccentric  and  the  alignment  of  the  E-cores  and  has  been  found  to  change  following 
reassembly  of  the  bearing.  A systematic  investigation  of  this  interaction  effect  was 
conducted  to  provide  the  variation  of  the  measured  pitch  and  yaw  angle  with  roll  position. 
Some  typical  results  of  this  investigation  are  shown  in  Figs.  9 through  12.  Data  were 
taken  of  the  variation  in  i//  and  0 as  a function  of  <!>  and  the  total  angle  in  the  vertical 
plane  for  sting  roll  positions  in  15-dcg  increments.  Accuracy  of  the  measurements  is  within 
±0.01  deg  in  0 and  i/^.The  sting  roll  position,  and  total  angle  of  attack,  6,  are  related 
to  the  Euler  angles  for  0 = 0 as  follows: 


5 = \js\n~  lif  I-  c,os^  ip  sin^  6 (36) 


lanQ  = y (37) 

uris  lir  sin  0 
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Figs.  9 through  1 2 reveal  some  important  facts  concerning  the  roll  interaction  effect  upon 
6 and  ^ measurements.  With  reference  to  these  figures,  the  magnitude  of  the  correction 
approaches  0.1  deg  for  practically  all  combinations  of  S and  fi.  In  addition,  there  is 
significant  variation  of  the  correction  as  a function  of  ^ with  no  apparent  symmetry  about 
particular  sting  roll  positions.  If  the  interaction  was  caused  by  imperfections  in  the  eccentric 
only,  symmetry  should  have  appeared  as  evidenced  by  matching  corrections  as  a function 
of  0 in  90-degree  sting  roll  positions  (compare  Figs.  9,11,  and  1 2).  Apparently  the  E-core, 
bearing  alignment,  and  eccentric  imperfections  all  play  a major  part  in  the  lack  of  symmetry 
of  the  roll  interaction  effect. 


In  addition  to  the  modeling  or  process  noise  resulting  from  the  roll  interaction  effect 
and  from  the  methods  used  for  bearing  calibration,  process  noise  due  to  bearing  moments 
can  also  be  present  in  typical  3DOF  data.  These  bearing  moments  are  small  and  highly 
nonlinear  and  vary  with  the  bearing  gas  pressure.  A special  pendulum  designed  for  the 
gas  bearing  was  used  to  study  nonlinear  bearing  moments  and  their  effect  on  3DOF  motion. 
Assuming  that  no  aerodynamic  moments  act  on  the  pendulum,  the  only  moment  caused 
by  the  pendulum  is  that  resulting  from  the  eg  displacement  from  the  vertical  axis  as 
depicted  in  Fig.  13.  Then,  in  the  aero  ballistic  system,  it  may  be  shown  that 


= — W / cos  sin  9 


A 

M = -W  / sin  Ilf 
z 


(38) 


where  2 is  the  distance  between  the  bearing  pivot  point  and  the  pendulum  center  of  mass 
and  W is  the  pendulum  weight.  Comparison  of  Eq.  (38)  with  Eqs.  (9)  and  (16)  yields 
a]  = -W2/q^Sd.  Geometric  properties  of  the  pendulum  and  assumed  conditions  for  data 
analysis  are  listed  below. 


Weight  = 10.38  Ibf,  I^  = 0.0305  slugs-ft^ 


ly  = Iz  = 0.04747  slugs-ft2,  2 = 0.058  ft 
S = 0.3975  ft2,  d = 0.7092  ft 


U.  = 1 ft/sec,  q,  = 0.5  Ibf/ftZ 

Experimental  (p,  6,  and  ^ data  were  taken  using  the  configuration  as  shown  m Fig. 
1 3 for  variations  in  the  bearing  pressure,  total  angle  of  attack  of  the  pendulum,  and  roll 
rate.  For  low  bearing  pressure  equal  to  100  psi,  analysis  of  the  3DOF  data  resulted  in 
excellent  fit  of  the  linear  moment  distribution  as  given  in  Eq.  (38).  Typical  plots  of  the 
resulting  solution  and  the  residuals  in  both  6 and  \i/  are  shown  in  Fig.  14.  Analysis  was 
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performed  employing  the  symmetric  version  of  the  coefficient  extraction  code,  Table  1, 
yielding  an  equivalent  Cm  a of  -3.70  1/rad.  With  reference  to  Fig.  14,  all  residuals  are  of  the 
same  level  as  measurement  errors  caused  by  the  roll  interaction  effect  discussed  earlier. 
In  addition,  a small  rolling-moment  coefficient,  Cgo  = -0.79  x 10-^  1/rad,  was  extracted 
from  the  data. 

As  the  bearing  pressure  was  increased,  the  resulting  3DOF  motion  became  highly 
nonlinear.  Large  differences  in  the  equivalent  Cm  „ for  small  compared  to  moderate  angles 
of  attack  were  detected.  The  extracted  values  of  the  equivalent  Cm  a for  the  moderate 
and  small  angles  of  attack  were  -3.79  1 /rad  and  -1 .32  1/rad,  respectively.  The  corresponding 
roll  moment  coefficient  increased  from  -0.000079  1/rad  for  moderate  angles  to  -0.00121 
1/rad  for  small  angles.  Thus,  as  bearing  pressure  is  increased  from  100  to  400  psi,  the 
bearing-induced  roll  torque  increases  by  an  order  of  magnitude  from  -1.1  x 10'^  ft-lbf 
to  -1.7  X 10-4  ft.ibf. 

Figures  15  and  16  show  the  fit  of  the  computed  solution  to  the  experimental  data 
for  a bearing  pressure  of  400  psi  at  moderate  and  small  angles  of  attack,  respectively. 
The  moderate  angle-of-attack  data  exhibited  large  deviations  from  motion  produced  by 
the  moments  as  pven  in  Eq.  (38).  This  is  evident  when  the  magnitude  and  character 
of  the  residuals  as  given  in  Fig,  1 5 are  compared  to  those  for  low  bearing  pressure.  Fig. 
14.  The  results  of  the  small  angle-of-attack  data.  Fig.  16,  indicate  good  agreement  between 
the  data  and  the  computed  solution  based  on  the  moment  expressions  as  given  in  Eq. 
(38);  however,  the  magnitude  of  the  equivalent  Cm  a is  a factor  of  three  smaller  than  the 
expected  value. 

Thus,  for  normal  operation  of  the  gas  bearing  at  bearing  pressures  of  400  psi,  nonlinear 
moments,  near  yp  and  6 equal  to  zero,  induced  by  the  bearing  are  present  in  typical  3DOF 
data  with  magnitudes  proportional  to  My  = 0.330  a ft-lbf  and  M,.  = -0.330  ]3  ft-lbf.  As 
bearing  pressure  is  decreased  from  400  psi,  nonlinear  bearing-induced  moments  are 
substantially  reduced.  For  a Mach  number  4 test  environment  with  q^  = 460  psf  and 
a slightly  blunted  sphere-cone  with  reference  dimensions  of  d = 0.7  ft  and  S = 0.35  ft^, 
the  contribution  of  the  bearing-induced  moment  at  400  psi  operating  pressure  is  on  the 
order  of  1 percent  of  the  model  aerodynamic  pitching-moment  coefficient  with  Cma  = -0-22 
1/rad. 

The  previous  discussion  has  revealed  the  existence  of  pitch  and  yaw  measurement 
errors  caused  primarily  by  the  roll  interaction  effect.  In  addition,  there  are  errors  present 
in  the  roll  position  measurement.  The  accuracy  with  which  roll  position  can  be  measured 
is  directly  influenced  by  the  number  of  alternating  reflective  and  nonreflective  lines  on 
the  eccentric  (see  Figs.  7 and  8).  Currently,  a change  in  roll  orientation  is  sensed  every 
two  degrees  by  the  digital  system.  Thus,  for  moderately  high  roll  rates  compared  to  the 
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median  frequency  in  the  pitch-yaw  plane,  a roll  uncertainty  of  ± 1 deg  appears  in  the 
data  as  approximately  random  noise.  As  the  roll  rate  is  decreased  toward  the  median 
pitch-yaw  plane  frequencies,  the  roll  uncertainty  appears  less  random  and  can  introduce 
significant  errors  in  extracted  coefficients  depending  on  the  duration  of  the  data  and  the 
character  of  the  motion  to  be  fitted. 

5.3  TUNNEL  A EXPERIMENTAL  RESULTS 

In  Sections  5.1  and  5.2,  errors  associated  with  the  gas  bearing  measurement  system 
were  examined.  The  ability  to  extract  correct  aerodynamic  coefficients  (up  to  ±0.1  deg 
in  pitch  and  yaw)  in  the  presence  of  random  measurement  noise  caused  by  the  gas 
bearing  was  shown  to  exist.  Even  though  these  errors  are  not  random,  successful 
analysis  of  low  bearing  pressure  bench  test  data  indicates  that  these  measurement  errors 
for  linear  systems  may  be  treated  as  random  measurement  noise.  The  successful  application 
of  the  parameter  identification  method  as  developed  here,  however,  also  strongly  depends 
on  the  precise  definition  of  process  noise  associated  with  tunnel  flow  nonuniformities. 
Depending  upon  the  exact  form  of  these  flow  nonuniformities,  extracted  aerodynamic 
parameters  with  large  bias  may  result.  In  order  to  examine  the  ability  of  the  present 
method  to  extract  "correct"  aerodynamic  coefficients  from  typical  wind  tunnel  results, 
3DOF  data  were  taken  in  the  AEDC-V-KF  Tunnel  A at  a Mach  number  and  Reynolds 
number  of  4.02  and  3.6  million  per  foot,  respectively. 

Tunnel  A is  a continuous,  closed-circuit,  variable  density  wind  tunnel  with  an 
automatically  driven  flexible-plate-type  nozzle  and  a 40-  by  40-in.  test  section  (Fig.  17). 
The  tunnel  can  be  operated  at  Mach  numbers  from  1.5  to  6.0  at  stagnation  pressures 
from  29  to  200  psia,  respectively,  and  at  stagnation  temperatures  up  to  750°R 
depending  upon  Mach  number  and  pressure  level.  Minimum  operating  pressures  range 
from  about  one-tenth  to  one-twentieth  of  the  maximum  at  each  Mach  number.  The 
tunnel  is  equipped  with  a model  injection  system  which  allows  removal  of  the  model 
from  the  test  section  while  the  tunnel  remains  in  operation. 

The  configurations  tested  were  a 6.7-deg  sharp  and  6.7-deg  spherically  blunted  cone 
with  mass  and  geometric  properties  as  listed  below  (see  Fig,  18). 

d = 0.6667  ft,  S = 0.3491  ft2 

eg  location  = 12.00  in.  from  model  base 

sharp  cone:  Ix  = 0.0301  slugs-ft^,  ly  = 0,2583  slugs-ft^,  = 0,2565  slugs-ft^ 

sphere  cone:  Ix  = 0.0304  slugs-ft^,  ly  = 0.2519  slugs-ft^,  1^  = 0.2514  slugs-ft^, 
dn/d  = 0.2275 
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These  two  configurations  do  not  sufficiently  challenge  the  analysis  capabilities  of  the 
asymmetric  or  symmetric  program  in  terms  of  model  aerodynamics.  They  do,  however, 
provide  baseline  data  on  bodies  whose  aerodynamics  are  known.  This  information  may 
be  used  to  assess  the  influence  of  small  tunnel  flow-field  nonunifoimities  on  the  sensitivity 
of  the  coefficient  extraction  technique. 

The  data  analyzed  are  shown  in  Figs.  19  and  20  where  the  variables  0,  0,  and 
are  plotted  against  time.  Figure  19  shows  typical  data  for  the  sharp  cone  model. 
Sphere-cone  results  are  given  in  Fig.  20  with  initial  conditions  similar  to  those  for  the 
sharp  cone.  Comparison  of  0 and  ^ versus  time  in  Fig.  19  shows  unexpected  frequency 
variations  and  damping  effects  normally  not  associated  with  motion  of  a sharp  cone.  For 
the  flow-field  and  body  constraints  listed  above,  analysis  of  the  sharp  cone  data  should 
yield  symmetric  aerodynamic  coefficients  with 

~ Cm  a ” -0.23  1/rad,  CmagZ  ~ ■Cnp5  2 ~ -1.13  1/rad^ 

Cnr  = Cmq  = -5.5  1 /rad,  Cgp  = -0.0081  1/rad  (Ref.  20) 

Analysis  of  the  data  as  given  in  Fig.  19  was  made  using  the  symmetric  program  assuming 
linear  aerodynamics  and  no  flow-field-induced  process  noise.  Analysis  of  8.44  sec  of  the 
data  resulted  in  extracted  coefficients  with 

Cma  = -0.225  1/rad,  Cmq  = -6.0  1/rad 
with  RMS  deviations  in  the  coefficients  of 

E(Cma)  = 0.0003  1/rad,  E(Cmq)  = 0.23  1/rad 

Comparison  of  the  extracted  with  the  corresponding  expected  value  of  each  of  the 
coefficients  reveals  good  agreement  which  is  reflected  in  the  size  of  the  calculated  RMS 
deviations.  However,  examination  of  the  residuals  and  the  fit  of  the  calculated  solution 
to  the  experimental  data  as  shown  in  Fig.  2 1 does  not  substantiate  the  validity  of  these 
coefficients. 

Correct  aerodynamic  coefficients  will  not  be  extracted  from  the  sharp  cone  data  until 
the  process  noise  in  the  experimental  data  is  properly  modeled.  Tunnel-fixed  flow 
nonuniformitics  arc  the  most  probable  source  of  the  process  noise  causing  the  poor  fit 
of  calculated  and  experimental  data.  Thus,  as  a first  approximation  in  the  modeling  of 
the  process  noise,  tunnel-fixed  stiffness  and  damping  moments  that  may  differ  in  both 
the  X,  Y and  X,  Z planes  are  assumed  as  shown  in  Eq.  (39). 
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M ,,  = q S(l(t|0  + t„  0) 
y * ^oc  ^ 2 ij  * 

. (39) 

Mzt  = ^ ‘4^^ 

oo 

In  the  aeroballistic  and  body-fixed  axis  systems,  the  moments  as  given  in  Eq.  (39)  transform 
to 


A 

Mx 

= MyT 

sin  \J/  cos  0 - 

M sin  0 

A 

= MyT 

cos  i/r 

A 

Mz 

= My,. 

sin  ifj  sin  0 4- 

My,,  cos  0 

= 

MyT 

sin  ijj 

COS  0 - 

M^,|,  sin  0 

M,  = 

MyT 

cos  i/r 

cos  0 4 

(My,,  sin  i/r  sin  0 4-  M^j  cos  0)  sin  (j> 

M,  = 

-My 

rp  cos 

[f/  sin  (f) 

+ (Myj  sin  ^ 

sin  '0  4-  My,,  cos  0)  cos  0 

Thus,  tunnel-fixed  flow^  effects  as  described  by  Eq,  (39)  may  be  extracted  from  3DOF 
data  using  the  symmetric  or  asymmetric  programs  by  including  the  moment  expansions 
as  given  in  Eqs.  (40)  and  (41),  respectively.  In  the  presence  of  other  unknown  process 
noise,  it  is  not  possible  to  extract  the  coefficients  ti,  t2,  ts,  and  t4  and  the  linear 
aerodynamic  damping  and  pitching-moment  coefficients  simultaneously  from  3DOF  data. 
This  is  due  to  the  overspecification  of  the  linear  stiffness  and  damping  effects  in  the 
moment  expansions.  This  may  be  shown  by  considering  the  linear  stiffness  terms  for 
moments  about  the  y axis  in  the  body-axis  system.  A parallel  derivation  may  be  shown 
for  the  damping  terms.  If  one  substitutes  Eq.  (39)  into  Eq.  (41),  neglecting  terms  involving 
sin  il/  sin  d,  and  uses  the  results  of  Eq.  (14),  Eq.  (42)  results,  as  follows: 

“ '■1  ^ li  + t/f  cos  0 sin  6 (42) 

+ I (sin  ij/  sin  + cos  tj/  sin  0 cos  + ... 

If  one  expands  terms  involving  ^ and  0 in  Eq.  (42)  and  rearranges,  Eq.  (42)  reduces 
to 

= (t,  X a^j)  0 cos<^  .r  (tg  + a^j)  ^ sin  <;!.  + ...  ^43^ 

Examination  of  Eq.  (43)  shows  that  the  coefficients  am  1 , ti , and  ts  for  small  angular 
motion  are  not  independent.  Thus,  only  two  of  the  three  parameters  (am  1,  ti , and  ) 
can  be  extracted  from  3 DOF  data  simultaneously.  Attempts  at  extracting  all  three 
coefficients  will  lead  to  divergence  of  the  iteration  scheme  in  the  parameter  extraction 
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technique.  As  stated  previously,  similar  constraints  on  the  stiffness  and  damping  coefficients 
in  the  aeroballistic  axis  system  and  on  the  damping  coefficients  in  the  body-fixed  axis 
system  may  be  shown  to  apply. 

In  addition  to  the  tunnel-fixed  stiffness  and  damping  effects,  constant  bias  in  the 
measured  pitch  and  yaw  displacements  may  also  be  present.  These  can  be  caused  by 
imprecise  alignment  of  the  zero  yaw/pitch  axes  with  respect  to  the  tunnel-fixed  axis  and/or 
flow  nonuniformities  causing  tunnel-fixed  trim  moments.  This  effect  has  been  included 
in  the  form  of  artificial  flow  angularity  corrections  in  the  pitch  and  yaw  planes  (0x  and 
i//i , respectively). 

The  asymmetric  parameter  identification  program  presented  in  Table  2 was  employed 
to  extract  the  tunnel-fixed  process  noise  (as  described  abovej  from  the  sharp  cone  data 
(Fig.  19).  Results  of  the  analysis  are  shown  in  Table  7.  As  indicated  in  the  table,  only 
the  tunnel-fixed  process  noise,  Ij,x,  Cgo,  and  the  boundary  conditions  were  allowed  to 
vary.  Due  to  the  constraints  on  the  tunnel-fixed  process  noise,  the  model  aerodynamics 
were  fixed  at  theoretical  estimates  given  previously.  The  calculated  solution,  residuals,  and 
fit  of  the  computed  solution  to  the  experimental  data  are  shown  in  Figs.  22  and  23 
for  two  data  sample  periods. 

Analysis  of  the  results  presented  in  Table  7 indicates  the  existence  of  tunnel  stiffness 
effects  in  both  the  pitch  and  yaw  planes.  In  addition,  a small  product  of  inertia  was 
extracted  from  the  data.  Tunnel-fixed  damping  effects  were  also  extracted,  indicating 
decreasing  magnitudes  for  increasing  data  sample  intervals.  This  effect  may  be  caused  by 
higher  order  coefficients  in  the  tunnel-fixed  damping  terms.  Examination  of  Fig.  22  shows 
excellent  fit  of  the  computed  to  the  experimental  data  with  residuals  in  6 and  V'  bounded 
by  ±0.1  deg.  The  best  fit  of  the  experimental  data  to  the  calculated  solution  is  bounded 
by  the  measurement  errors  caused  by  the  gas  bearing  as  discussed  in  Section  5.2.  Doubling 
of  the  data  sample  interval  as  shown  in  Fig.  23  resulted  in  larger  residuals  with  a dominate 
frequency  approximately  the  same  as  the  fundamental  frequency  of  the  experimental  data 
(compare  Figs.  23c  and  d with  Figs.  19a  and  b,  respectively).  This  is  most  probably  an 
effect  of  frequency  mismatch  caused  by  higher  order  terms  not  included  in  the  tunnel 
stiffness  model,  Eq.  (39).  Considering  the  results  given  in  Table  7 and  depicted  in  Figs. 
22  and  23,  the  fit  of  the  calculated  3DOF  solution  to  the  experimental  sharp  cone  data 
is  significantly  improved  over  that  given  in  Fig.  21.  Thus,  tunnel-fixed  process  noise 
governed  by  Eq.  (39)  and  Izx  appear  to  be  valid  contributions  to  the  parameters  to  be 
evaluated. 

In  order  to  further  evaluate  the  tunnel-fixed  process  noise  and  the  ability  of  the 
asymmetric  program  to  extract  model  aerodynamics  in  the  presence  of  this  process  noise, 
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the  sphere  cone  data  as  shown  in  Fig.  20  were  analyzed.  Using  the  asymmetric  program, 
aerodynamic  moment  coefficients  and  tunnel-fixed  process  noise  were  extracted  from  the 
data  as  given  in  Fig.  20  and  are  shown  in  Table  8.  Since  model  aerodynamic  coefficients 
were  allowed  to  vary,  one  damping  and  one  stiffness  coefficient  describing  the  tunnel-fixed 
effects  were  held  constant.  Considering  the  RMS  deviations  and  values  of  tj,  t2,  ta,  and 
t4  as  given  in  Table  7,  the  stiffness  term,  ti  and  damping  term,  t2,  were  assumed  constant 
and  equal  to  0.017  1/rad  and  0,  respectively.  Examination  of  Table  8 shows  large  variations 
in  the  extracted  aerodynamic  coefficients  for  different  data  sample  intervals.  Moreover, 
since  the  sphere  cone  is  axially  symmetric,  the  relationship  given  in  Eq.  (16)  should  hold. 
As  the  data  sample  interval  increases  in  Table  8,  the  extracted  coefficients  tend  to  become 
less  like  those  caused  by  an  axisymmetric  body.  These  effects  are  further  illustrated  by 
comparison  of  the  graphical  results  in  Figs.  24,  25,  and  26  which  correspond  to  this  table. 
Thus,  in  addition  to  those  sources  of  process  noise  described  here,  others  must  exist  which 
cause  the  significant  disagreement  in  the  generated  solution  and  experimental  data.  Even 
though  the  preceding  analysis  of  the  sphere  cone  data  as  given  in  Fig.  20  does  not  yield 
acceptable  extracted  aerodynamic  coefficients,  it  does  indicate  some  problem  areas 
encountered  in  the  analysis  of  3DOF  data  taken  in  continuous  wind  tunnels  as  well  as 
concepts  for  establishing  process  noise  models. 

6.0  CONCLUDING  REMARKS 

The  previous  sections  have  presented  the  mathematical  development  of  coefficient 
extraction  programs  and  practical  application  of  the  programs  in  continuous  wind  tunnel 
flow  fields.  The  programs  developed  here  are  based  on  the  principals  governing  output 
error  methods.  As  such,  extracted  parameters  are  highly  sensitive  to  measurement  and 
process  noise.  Thus,  identification  of  sources  of  both  process  and  measurement  noise  and 
modeling  of  these  effects  were  extensively  examined. 

Noise  attributable  to  the  gas  bearing,  tunnel  flow  nonuniformities,  and  inappropriate 
aerodynamic  modeling  was  catagorized  as  being  the  most  important.  Measurement  noise 
caused  by  the  gas  bearing  was  shown  to  contribute  errors  in  i//  and  0 measurements  up 
to  ±0.1  deg  with  the  magnitude  of  these  errors  dependent  upon  the  value  of  \/),  0,  and 
(p  and  the  gas  bearing  pressure.  In  addition,  small,  bearing-induced  moments  were  also 
identified.  Even  though  the  magnitude  of  the  bearing-induced  errors  can  be  large,  it  was 
concluded  that  these  errors  for  most  situations  could  be  modeled  as  random  noise.  This 
allows  for  accurate  parameter  identification,  assuming  sufficiently  large  data  samples  and 
no  other  noise  sources. 

In  order  to  define  the  effects  of  process  noise  on  the  extracted  coefficients,  sharp 
and  spherically  blunted  cone  data  taken  in  Tunnel  A were  analyzed. 
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It  was  shown  that  the  process  noise  inherent  in  the  sharp  cone  and  sphere  cone 
data  was  large  and  could  not  be  treated  satisfactorily  as  linear,  tunnel-fixed  stiffness  and 
damping  moments  independent  of  the  model  configuration.  In  addition,  this  process  noise 
was  of  sufficient  level  to  cause  significant  error  in  the  extracted  coefficients.  The  analysis 
of  the  sharp  cone  and  sphere  cone  data  did  provide,  however,  concepts  for  establishing 
process  noise  models  in  similar  situations.  Further  studies  are  planned  on  the  definition 
of  process  noise  in  Tunnel  A and  on  the  effects  of  these  noise  sources  on  extracted 
aerodynamic  coefficients.  Also,  more  work  is  planned  on  the  application  methods  of  the 
programs  developed  herein. 
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Body-Fixed 
Coordinate  System 


Figure  1.  Arbitrary  body  showing  axis  system  for  3DOF 
equation  development. 
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Figure  2.  Plane-fixed  coordinate  systen^  showing  Euler  angles  and 
velocity  components  relative  to  space-fixed  system. 
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a.  Symmetry  operation  for  mirror  symmetry 
in  X,  z plane 


/ 

X 


b.  Symmetry  operation  for  90-deg  rotational 
symmetry  about  the  x axis 

Figure  3.  Physical  description  of  symmetry  operations  on 
axis  systems. 
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Figure  4.  Flow  Chart:  symmetric  and  asymmetric  parameter  extraction  program. 
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d.  Residual  in  \j/  versus  time 
Figure  5.  Continued. 
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Calculated 


f.  Fit  of  calculated  solution  to 
experimental  data  in  \l/,  6 
.plane,  time  increment  1.98 
to  3.98  sec 


g.  Fit  of  calculated  solution  to 
experimental  data  in 
plane,  time  increment  4.98 
to  5.98  sec 
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Figure  5.  Concluded. 
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Figure  8.  3DOF  balance  schematic, 
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a.  6 correction  versus  0 

iction  effect  for  sting  roll  position  £2  = 0 deg. 
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b.  rp  correction  versus  0 
Figure  11.  Concluded. 
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b.  Computed  }p  versus  time 
Figure  14.  Continued. 
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d.  Residual  in  *//  versus  time 
Figure  14.  Continued. 
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g.  Fit  of  calculated  solution  to  experimental  data 
in  9 plane,  time  increment  6.36  to  8.6  sec 
Figure  14.  Concluded. 
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b.  Computed  ^ versus  time 
Figure  15.  Continued. 
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d.  Residual  in  ^ versus  time 
Figure  15.  Continued. 
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a.  Computed  d versus  time 


b.  Computed  <//  versus  time 

Figure  16.  Bench  test  results  for  small  angle  of  attack  and  400-psi  bearing  pressure. 
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b.  i|/  versus  time 
Figure  19.  Continued. 
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Figure  21.  Graphical  solution  of  sharp  cone  data  given  in  Fig.  19  for  a fit 
interval  of  8.44  sec  with  no  tunnel-fixed  process  noise. 


H -I 

11.00  12-00 


AEDC-TR -78-10 


AEDC-TR-78-10 


AEDC-TR-78-10 


e.  Fit  of  calculated  solution  to  experimental  data  f.  Fi1 

in  <//,  0 plane,  time  increment  0 to  0.98  sec  in 

Figure  21.  Continued. 


f.  Fit  of  calculated  solution  to  experimental  data 
in  d plane,  time  increment  0.98  to  1.98  sec 
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k.  Fit  of  calculated  solution  to  experimental  data  I.  Fit  of  calculated  solution  to  experimental  data 

in  \p,  6 plane,  time  increment  5.98  to  6.98  sec  in  \j/,  6 plane,  time  increment  6.98  to  7.98  sec 

Figure  21.  Concluded. 
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Figure  22.  Graphical  solution  of  sharp  cone  data  given  in  Fig.  19  for  a 
fit  interval  of  4.4  sec  with  fixed  aerodynamics. 
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d.  Residual  in  \p  versus  time 
Figure  23.  Continued. 
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k.  Fit  of  calculated  solution  to  experimental  data  I.  Fit  of  calculated  solution  to  experimental  data 

in  yji.  6 plane,  time  increment  558  to  6.98  sec  in  \j/,  6 plane,  time  increment  658  to  7.98  sec 

Figure  23.  Concluded. 
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d.  Residual  in  ^ versus  time 
Figure  24.  Continued. 
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e.  Fit  of  calculated  solution  to  experimental  data  f.  Fit  of  calculated  solution  to  experimental  data 

in  0 plane,  time  increment  0 to  0.68  sec  in  G plane,  time  increment  0.68  to  1.38  sec 


Figure  24.  Continued. 
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j.  Fit  of  calculated  solution  to  experimental  data 
in  \1/.  0 plane,  time  increment  3.48  to  4.38  sec 


Figure  24.  Concluded. 


k.  Fit  of  calculated  solution  to  experimental  data 
in  \Ij,  0 plane,  time  increment  4.38  to  5.16  sec 
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a.  Computed  6 versus  time 

Figure  25.  Graphical  solution  of  sphere  cone  data  given  in  Fig.  20  for  a 
fit  interval  of  8.32  sec. 
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c.  Residual  in  6 versus  time 
Figure  25.  Continued. 
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g.  Fit  of  calculated  solution  to  experimental  data  h.  Fit  of  calculated  solution  to  experimental  data 

in  <//,  0 plane,  time  increment  1.38  to  2.08  sec  in  ip,  0 plane,  time  increment  2.08  to  2.78  sec 

Figure  25.  Continued. 
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Fit  of  calculated  solution  to  experimental  data 
in  i//,  0 plane,  time  increment  2.78  to  3.48  sec 

Figure  25. 


Figure  25.  Continued. 
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m.  Fit  of  calculated  solution  to  experimental  data  n.  Fit  of  calculated  solution  to  experimental  data 

in  \l/.  6 plane,  time  increment  5.58  to  62B  sec  in  i//,  d plane,  time  increment  628  to  6.98  sec 


Figure  25.  Continued. 


o.  Fit  of  calculated  solution  to  experimental  data  p.  Fit  of  calculated  solution  to  experimental  data 

in  \p,  6 plane,  time  increment  6.98  to  7.68  sec  in  <//,  6 plane,  time  increment  7.68  to  8.28  sec 

Figure  25.  Concluded. 
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a.  Computed  6 versus  time 

Figure  26.  Graphical  solution  of  sphere  cone  data  given  in  Fig.  20  for  a 
fit  interval  of  11.98  sec. 
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c.  Residual  in  0 versus  time 
Figure  26.  Continued. 
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d.  Residual  in  \]/  versus  time 
Figure  26.  Continued. 
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e.  Fit  of  calculated  solution  to  experimental  data  f.  Fit  of  calculated  solution  to  experimental  data 

in  yjj.d  plane,  time  increment  0 to  0.68  sec  in  t//.  0 plane,  time  increment  0.68  to  1.38  sec 

Figure  26.  Continued. 


g.  Fit  of  calculated  solution  to  experimental  data 
in  \l/.  0 plane,  time  increment  1.38  to  3.48  sec 


h.  Fit  of  calculated  solution  to  experimental  data 
in  6 plane,  time  increment  3.48  to  4.18  sec 


Figure  26.  Continued. 
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i.  Fit  of  calculated  solution  to  experimental  data  j-  Fit  of  calculated  solution  to  experimental  data 

in  \Ij,  e plane,  time  increment  4.18  to  6.28  sec  in  i//,  d plane,  time  increment  6.28  to  6.98  sec 

Figure  26.  Continued. 


k.  Fit  of  calculated  solution  to  experimental  data  I.  Fit  of  calculated  solution  to  experimental  data 

in  \}/.  e plane,  time  increment  6.98  to  7.98  sec  in  d plane,  time  increment  7.98  to  8.68  sec 

Figure  26.  Continued. 
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m.  Fit  of  calculated  solution  to  experimental  data  n.  Fit  of  calculated  solution  to  experimental  data 

in  \Ij,  6 plane,  time  increment  8.68  to  9.98  sec  in  <//,  0 plane,  time  increment  9.98  to  10.68  sec 

Figure  26.  Continued. 


o.  Fit  of  calculated  solution  to  experimental  data  p.  Fit  of  calculated  solution  to  experimental  data 

in  d plane,  time  increment  10.68  to  11.28  sec  in  6 plane,  time  increment  11.28  to  11.98  sec 

Figure  26.  Concluded. 
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Table  1.  Fortran  Listing:  Symmetric  Parameter  Extraction  Program 


ANGLES 

VARNER  VERSION  1 FEB.  1977 

COMPUTES  AERO  COEFFICIENTS  FOR  AXIALLY  SYMMETRIC  BODY 
AEROBALLISTIC  AXIS  SYSTEMt AERODYNAMIC  SYMMETRY 
EXACT  3 DOF  EQUATIONS  UTILIZED 


IMPLICIT  REALMS (A-HtO*Z> 

REALMS  IXtlY*MPRM,IPRM 
REALMS  ALISTOO) 

REAL*8  CARD (10) 

REAL*4  TS(IOOO) (PHIS (1000) (THTAS(IOOO) iPSIS(lOOO) 

DIMENSION  AUX  (6(186)  (Y(186)  (DERY(186)  (WXX  (30)  ( ACONOO)  •CT( 30(31) 
1(WX(30) (AACON(30) (CC (30(60) 
external  FCTfORKGS 

COMHON/MEAS/PHIM(1000) (THTAH(IOOO) (PSIM(IOOO) 

COMMON/SOL/TdOOO)  iPHKlOOO)  (THTA  ( 1000)  (PSI  ( 1000) 
COMMON/COEF/CON(30) (CCON(30) (NC(NCON(30) 

COMMON/MATRIX/C (30(60) tRESSUM(S16PHl(SlGTHA(SlGPSI 
COMMON/CONST/RAT 10(MPRM( 1PRM( A( lANGLE 
equivalence  ( AUX ( 1 ) (CT  < 1 ) I ( (CC  < 1 ) (C ( 1 ) ) ( 

1 (PHISI1)(PH1(1))((THTAS(1)(THTA(1))((PSIS(I)(PS1(1)) 

DATA  ALIST/*CMA  '('CNQ  •(•CHAD  *('CMPA  •(•CMA0*«2  •( 

»«CMQD*«2  •(*CMOB**2  *(«CHRAB  *('CMQ0  '('CLO  *('CLP0  *( 

*'CMPAD**2* ( ’CLP  *(»CLP0**2  •(•DUMMY  •(•DUMMY  •(•DUMMY  •( 

••dummy  ((tCNR  TUN  •(•CMO  TUN  •(•CMA  TUN  •(•CNB  TUN  •(•FL0AN6  T*( 

••FLOANG  P*(»PH1  (0)  •(•PHIOOT  •(•THTA(O)  •(•THTADOT  •(•PSKO)  •( 

••PSIDOT  •/ 

NAMELIST/NAM/NSNTCH(NT0T(NC(NITER(DTMIN(B0UN0(AREA(DIA(IX(IY(R0(U( 

1AAC0N.NC0N(NSESNH*SPHI(STHTA(SPSI(AMPHI(AMTHTA(AMPSI( 

2SAMPLEiNDELTA(IUP0AT(INLC0N(IANGLE(DTGEN 

3(I6PnO(18LKNO(1RC6EN(SIGPHI(SI6THA(SIGPS1(NSTART(IRCSTR 

100  FDRMAT(1H1(20X(^»»**«***  INPUT  DATA  *•*••*•••(/( 

11X(*GR0UP  number  ■•(I5(SX(*BL0CK  number  ••iI5(/( 

19X('TIME*(10X('PSP(10X(«THETA^(10X('PHP(//) 

101  FORMAT(7F10.0) 

102  F0RMAT(5X(A(E13(6(1X) ) 

103  F0RHAT(5X(A8.2X(2(E13.6(SX)) 

lOA  FORMAT(/(1X(^1TERATION  N0<(I5(/(19X(*C0N^(13X(tDELTA  CON*) 

105  FORMAT(///(10X(*FINAL  CONSTANTS  AND  PROBABLE  ERRORS (10X(/) 

106  F0RMAT(5X(^C('(13(*)b^(E13.6(5X(^ ( • (E13.6( • ) • (2X( A8) 

10?  F0RMAT(1M1(5X(^0ATA  fit  COMPARISON^(//(11X(  »TIME't 

13X(*PSI/F1T«(3X(^PSI/EXP^(3X('0EL  PSI • .2X( •THTA/FIT^ (2X( •THTA/EXP* 
2(2X(^DEL  TMTA*(2X(*PHl/FIT*(3X(^PHI/EXP^(3X(^DEL  PHP(/) 

108  F0RMAT(SX(F10.6(9F10.3) 

109  format  I1X(****«  convergence  criteria  not  SATISFIED  AFTER*(IA( 
11X('ITERATI0NS  *••••) 

110  F0RMAT(1H1(1X(*SUM  OF  RESIDUALS* (E13.6(SX( •NEW  PROBABLE  ERROR  «•( 
1E13(6(5X('0LD  probable  ERROR  ••(E13.6) 

111  F0RMAT(7F10.4) 

112  FORMAT(1X(70(^*^)(/(1X(^******»*  ERROR  •**••••*• (/(1X( 

PNO  OF  CONSTANTS  ALLOWED  TO  VARY  EXCEEDS  TOTAL  NO  OF  DATA  POINTS 
2USEO^(/(1X(70(^*')) 
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Table  1.  Continued 


200  FORHAT(1Xi30('*') *1  AERODYNAMIC  COEFFICIENT  EXTRACTION  BY  CHAPMAN 

IKIRK  method  •t30(«*')«///«lXt'N5WTCH  sQ*  NORMAL  DATA  REOU 

2CTlON*«/tl7Xt'«lf  3D0F  DATA  GENERATEDt  NO  ITERATION' t/t 
317x*f>2*  3D0F  DATA  GENERATED  AND  STORED  ON  7S  PACK» • f/tlX« • lUPDAT 
Aa'.lS.KsOt  INITIAL  GUESSES  NOT  UPDATED  FOR  FOLLOWING  SHOT* t/t 
S17Xt'NE  Ot  INITIAL  GUESSES  UPDATED)*) 

201  FORMAT(//tlXt ’PARAMETERS  *•••*♦•* t/tlXt 'NO  OF  POINTS  IN)  "'tlSt 

1 'TOTAL  NO  VARIABLES  (NO  "' t IStSXt 'TOTAL  NO  ITERATIONS  (NITER) >'t 
215i/tlXf 'INTEGRATION  STEP  SIZE (DT) s' tE8.2t5Xt 'CONVERGENCE  BOUND  ON 
3 iteration  (B0UN0)B'tE8.2t//tlXt*B00Y  •••••««•' t/tlXt 'AREA  >'t 
♦Fl0.6t'Fr*»2't5Xt'OlA  s' tF10,6t 'FT* t/t IXt • IX  »* tFl2.7t *SLUGS-FT**2 
S'tSXt 'lY  s'tF12.7t'SLUGS-FT**2*) 

202  FORMAT(/tlXt'FR£E  STREAM  ♦***••**' t/tlXt 'DENSITY  (RO)  ■*tEHt7t 
1'SLUGS/FT**3* tSXt'U  a ' tFlO. A t'FT/SEC't//tlXt' COEFFICIENTS  ♦*•*•*•) 

203  format  (lAt'ACONCtlSt*  I aitE13.6t'  I • t A8t ' ) • tSXt 'NCON ( * t I3t ' ) **  1 12) 

301  FORMAT(lHl) 

302  FORMATdXt'lANGLE  at«I5t'(>0t  ALPA  AND  BETA  IN  TERMS  OF  VELOCITY  R 
lATlOSt't/tlTXt'NE  Ot  ALPA  AND  BETA  IN  TERMS  OF  ANGLE  OF  ATTACK  AND 

3 SIDESLIP) 't/tlXt'INLCON  a',l5t'(a  Ot  INITIAL  CONDITIONS  EVALUATED 

4 FROM  0ATA't/tI7Xt'NE  OtlNPUT  VALUES  USED)') 

303  FORMATdXt 'BEGINNING  RECORD  NUMBER  a',l5) 

401  FORMATdOAS) 

402  FORMATdXtlOAB) 

404  FORMATdXt  INPUT  CARO  IMAGES  ••<"t//) 

406  format (IHl) 

407  FORMATdXt*  NOTE  • 'tlOAB) 

502  F0RMATdXt69("")  t/tlXt'****  300F  GENERATED  DATA  STORED  BEGINNING 
20N  RECORD  NUMBER* t I5t2Xt '••••' t/t 1 Xt 'UNDER  IGPNO  a',I5t2Xt 
3*  AND  IBLKNO  >' t I5t /t IX t • I lal ,T * t /t IX t • I I >2tPHI < t/ t IXt 
4*Ila3,THrA't/tlXt*IIa4tPSI*t/tlXt69("") ) 

WRlTE(06t404) 

405  READ(05t40ltENDB403)CARD 
WRlTE(06t402)CARO 

60  TO  405 

403  REWIND  S 
READ(05t401)CARD 
wRlTE(06t406) 

DEFINE  file  07(1140t4016tLtIREC) 

INITIALIZATION 

NPRlNTaO 

IPRiNTaO 

998  REAO(05tNAMtENOB999) 

NaNTOT 

WRlTE(06t301) 

WRITE (06t200) NSwTCHt IUPDAT 
WRITE (06t302> lANGLEt INLCON 
WRITE(06t407)CARO 

WRlTE(06t201)NTOTtNCtNITERtOTMINtBOUNOtAREAtDIAtlXtlY 

WRITE(06t202)ROtU 

DO  204  laltNC 

ACONd)aAACON(I) 

204  WRITE (06t203) 1 t ACON( I ) t ALIST (I ) t I tNCON ( I ) 

PI  « 3.1415920*0 
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A • PI/lb0t0»0 
IF(NSWTCM>8f8t7 


data  input 


7 continue 
T(l)»0.0*0 
00  10  I»2tNT0T 
T<I)»T(I-l)*DTeEN 
PHIM(I)«0,D40 
THTAM(I)aO.O*0 
10  PSlM(I)a0,O«0 
PHIM(l)s0.0-»0 
THTAM<1>  >0«D*0 
PSIM(U»0,0»0 
QO  TO  9 

6 continue 

IRECal 

500  HEAO(07*IREC)IGPT,1BUKiII*TS 
IF(I6PN0.NE.I6PT>  GO  TO  500 
IFdBLK.NE.IBLKNO)  GO  TO  500 
1F<11.NE.I>  60  TO  500 
HEAD ( 07 • IREC) IGPNOf IBLKNO* I 1 tPHIS 
READ (07* IREC ) IGPNOt IBLKNOt 1 I .THTAS 
READ (07* IREC ) IGPNO* IBLKNO* 1 1 tPSIS 


C 

C 

C 


data  alteration 


DO  50A  I«1*NT0T 
T(I)>TS(I) 

PH1H(I)«PHIS(1) 

THTAH(I>aTHTAS(l> 

50A  psiM(n«Psis(n 

IF (NPRlNT.EQ.16PN0tANDtlPRlNT.E0. IBLKNO) 
t.RlTE(06tlOO)16PNOiIBLKNO 


60  TO  A 


IRCTMPalREC-A 

WR1TE(06»303)  IRCTMP 

LOOP«0 

NLOOPal 

DO  A laltNTOT 

WR1TE(06*102)T(1) *PSIM(1) *THTAM( I) •PHIMIl) 
LOOPaLOOP*! 

KLOOPaNLOOP«50 
IF(LOOP.LT.KLOOP)  60  TO  A 
NLOOPaNLUOP*! 

WRITE (06* too) I6PN0.IBLKN0 


A continue 
NPRlNTalGPNO 
IPRiNTalBLKNO 

CALL  part (SAMPLE*NOELTA»N»NTOT#NSTART) 

IF  (NSESNH.EO.DCALL  NOISE  (SPHl.STHTAiSPSI 
lNtSAMPLE«NOELTA*NTOT) 

IF(INLCON.EQ.O)  call  ICS(ACON) 

DO  A6  laliN 
PSlM(l)  a PSIM(I)*A 


f AMPHl * AMTHTA* AHPSl * 
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THTAM(I)  ■ THTAM(n*A 
46  PHIM(I)  a PHIM(I)«A 
9 continue 
NSUMaO 

00  12  laltNC 
12  CCON(n>0.0*0 
DO  11  laltNC 

1F(NC0N(1) .EQ.1>  CCON(I)«1.0«0 
II  IF(NC0N(I) .EQ.ll  NSUMaNSUM*! 

SUM  a NSUM 

DEFINE  PARAMETERS  AND  TRANSFORM  COEFFICIENTS  IN  SUB  REG 

OINF  a RO«IU**2)/2«D«0 

RATIO  a IX/IY 

IPRM  a IY/(QINF*AREAaOIA) 

MPRM  a IPRH*2.D»0*U/DIA 
call  REGIACON) 

XNaN 

•-HECKa3.0«0aXN«SUM 
IF(CHECK.LE.0.O«0)  GO  TO  997 
NTaNC*! 

EaQ.D^O 

NRaO 

integration 

20  continue 

CALL  INT6RL(FCTtORKGStNtERRORiAUXtYtOERY«DTMIN) 
IF(NSWTCM,GT.0)  go  to  47 
IFCNSUM.EO.O)  GO  TO  47 

CHECK  FOR  CONVERGENCE 

NRaNR^l 

ElaUS0RT(RESSUM/(3«0«0*XN>SUM) ) 

RESSUMaOSQRT (RESSUM) 

WRITE (06tl 10) RESSUM*E1«E 
IF (DABS (El-E) -BOUND) 36*36*21 

21  IF(NR.NITER)22«22*36 

22  NTRlGaO 
38  CONTINUE 

IF(NSUH.EQ.l)  GO  TO  33 
E-El 

contract  C matrix  • REMOVE  ZERO  ROWS  AND  COLUMNS 

DO  23  Ial*NC 
DO  23  J-ltNT 

23  CT(1.J)-C(I*J) 

LLaO 

DO  90  lal.NC 

92  1F( (I*LL) .GT.NC)  GO  TO  90 
IF (NC0N(1*LL) >91*91*93 
91  LLbLL*1 
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60  TO  92 

93  L-UL*1 

00  94  K«lfNT 

94  C(l«K)aCT(LfKI 
90  CONTINUE 

LL«0 
LMaNSUM 
DO  SO  I«ltNC 

82  1F( (I«LL) .GT.NC)  GO  TO  SO 
1F(NC0N(I*LL} >81t81f83 

81  LL>LL*1 
60  TO  82 

83  LaLL^I 

DO  84  K>1*LM 
64  CT<KiI)«C(KfLI 

80  continue 

00  8S  laltLM 
85  CT(IfLM«l)aC(I*NT) 

LaNSUM 
LLaNSUM^l 
DO  71  Ial,L 
DO  71  jaliL 
71  C(ItJ)aCnitJ) 

MATKIX  inversion  OF  C 

lF(NTRlG)40i27t40 

27  CALL  INV<CTiLtLL«WXtCCI 
60  TO  41 

40  call  INV(CTtL*LtWX«CCI 

41  00  32  I>lfNC 
C(I*LL)«0.0«0 
WX(l)aCT (Xtl) 

32  MXX(l)a0»O*O 
LMaO 

expand  C inverse  to  include  zero  ROMS  AND  COLUMNS 

DO  31  laitNC 
1F(NC0N(I) I28t20t30 

28  LHaLM-1 
GO  TO  31 

30  lTaLM*I 
C(1*LLI*CT(IT«LL) 

MXX(l)aMX(IT) 

31  continue 
60  TO  34 

SOLVE  FOR  CORRECTION  TO  COEFFICIENTS 

33  DO  35  I«ltNC 
1F(NC0N(1)-I)39t49»49 

49  C(ItNT>aC(I«NT)/C(ItI) 

MXXlI>al.D*0/C(ItI) 

GO  TO  35 
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39  C<ltNT>>U.O«0 
WXX(1)«0.0*0 

35  continue 
LL«NT 

34  lF{NTRI6.EQ.n  60  TO  42 
WRITE(06tl04|NR 
DO  60  I«lfNC 
60  WXX(I)«C(ItLL) 

transform  COEFFICIENTS  AND  CORRECTIONS  IN  SUB  REGOUT 

CALL  REGOUTIACONfWXX) 

DO  50  I>lfNC 

SO  WRlTE(06tl03>ALIST(I) tACON(I) tWXXII) 

DO  43  I«1«NC 

43  CON(I)«CON(I)*C(I»LL) 

60  TO  20 

36  continue 
NTRlGsl 

37  GO  TO  38 
42  continue 

COMPUTE  PROBABLE  ERROR  IN  COEFFICIENTS 
DO  44  I>ltNC 

44  WXX<II>E1*0SQRT(DABS(WXX(I>)) 

CALC  REGOUT (ACONtMXX) 

WRlTE(06tl05) 

IF <NR. GT. NITER) WRITE (06i 109) nr 
DO  4S  I^ItNC 

45  WRlTE(06tl06)  I»AC0N(I) tWXX(I) .ALISTCI) 

47  WRITE(06*107) 

LOOP»0 

NLOOPbI 

IRECbIRCGEN 

0TTEMPsT(2)-T«n 

DO  48  I>1.N 

TT»T(I> 

Z1  • PSKD/A 
Z2  ■ PSIM(1)/A 
Z3  ■ THTA(I)/A 
Z4  • THTAM(I)/A 
Z5  • PHI (I) /A 
Z6  • PHIM(I)/A 
Z10-Z1-Z2 
Z11«Z3-Z4 
Z12«ZS«Z6 

WRITEI06»108)TT»ZltZ2tZ10tZ3«Z4tZlIfZ5tZ6»Z12 

THTASI1>«Z3 

PSIS<I)>Z1 

TS(I)»TT 

PHIS<1)«25 

IF(NSWTCH,NE.O)  GO  TO  597 

TS(1)>Z10 

PHIS(I)-Z11 
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597  continue 
LOOP»LOOP*1 
KLOOP»NLOOP«50 
1F(L00P«LT*KL00P>  60  TO  46 
NL00P>NL00P*1 
WR1TE(06»107) 

48  continue 

IF(NSWTCH.NE.O)  GO  TO  599 
C 

C COMPARISON  OF  FIT  DATA  AND  EXPERIMENT  STORED 
C 

IRECbIRCSTR 

lF(lRCSTH.GE,IRCeEN-3.AND.lRCSTR.LE.lRC6EN43)  GO  TO  599 
WRITE ( 06  »598»  NSWTCHi IRCSTRt IGPNO* laLKNOiOTTEMP 

598  FORHATIlHltlXf 'NSWTCH  »' tI5i/f IXi •3D0F  DATA  STORED  BEGINNING  AT  RE 
ICORD  i«15*2X« 'UNDER  I6PN0  ■< t IStZXt ' AND  IBLKNO  *I5f/»iXi • II»5|TH 

2TA»*/ilXf»Ila6tTHTA-THTAM't/tlXf*lI»7tPSI'*/ilX**II»8*PSI-PSIM'» 

3/flXt'DT  ■*fF7.4*2X«'SEC'tlHI) 

1I>5 

WR1TE(07*IRECI IGPNO.lBLKNOtlliTHTAS 
Il»6 

WRITE  (07 'IREOlGPNOi  IBLKNO*  I ItPHIS 
Il«7 

WRITE  (07 'IREOIGPNO*  IBLKNO*  1 1*PSIS 
11*8 

wRITE(07'1REC)1GPNO*IBLKNO*II*TS 

599  CONTINUE 

IF(NSWTCH.NE.2)  60  TO  501 
C 

C 300F  GENERATED  SOLUTION  STORED 
C 

WRITE(06*502) IRC6EN • I GPNO  * I BLKNO 
1I>1 

WRITE (07 'IRECI I6PNO* IBLKN0*II*TS 
Il«2 

WRITE(07*IRCC)IGPNO»IBLKNO*II*PHIS 

11-3 

WRITE ( 07 • IREC) 16PN0* IBLKNO* 1 1 *THTAS 
11-4 

WRITE ( 07 • IREC) IGPNO* IBLKNO* 1 1 *PSIS 
501  CONTINUE 
GO  TO  996 
997  WRITE(06*112) 

996  continue 

IF(IUPDAT.EQ.O)  60  TO  998 
DO  995  lal*NC 
995  AAC0N(I)«AC0N(1) 

60  TO  998 
999  continue 
STOP 
END 


SUBROUTINE  0RK6S(XI*Y*0ERY*N01M*FCT«AUX*DT*NDT) 
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c 

c 

c 


IMPLICIT  REAL*e<A>H»0<-Z) 
dimension  Y(l)tOERYIl)tAUX{6fl) 

XBXl 

DO  S JaltNOT 
DO  1 IslfNOIM 

1 AUXl6tnaY(I) 

CALL  FCTIXtYtOERY) 

DO  2 lalfNOIM 
AUX(l«I)aDERY(l) 

2 Y ( I ) >AUX  <6 f 1 ) * .5D»0*DT*AUX ( 1 • 1 1 
XTEMPaX*t50*0«DT 

CALL  FCT(XTEMP»YfDERY) 

DO  3 laltNOIM 
AUX(2tI)"DERY(I) 

J Y(I>BAUX(6f I)*.5D«0*DT*AUX<2f n 
call  FCTIXTEMP*Y*DERY) 

DO  A laltNOlH 
AUX(3«I)>DERY(I) 

A Y(I>aAUX(6iI)«0T*AUX(3tI) 

X«X*DT 

CALL  FCT(XtY*OERYI 
DO  5 lalfNOIM 
AUXU,l)aOERY(I) 

5 Y(I)aAUX(6tI>*0T*IAUX(ltl)*(AUX(2tI)*AUX(3tl) )*2tD*0 
1«AUX(A«I) )/6.0«0 
RETURN 
END 


SUBROUTINE  FCT (TtXiDERX) 

IMPLICIT  REALMS (A-H»0«Z> 

REAL*8  MX0lXtMY0IfMZ0I«MltM2«M3*MA6tMPRMfIPRM 
dimension  X(l)tOERX(l) 

COMMON/COEF/CON(30) tCCON(30)fNC*NCON(30) 

COMMON/CUNST/RATIO«MPRM«lPRMtAtK 

PSlaXCS) 

THTAaX(3> 

PHlax(l) 

PSIPaX<6) 

THTAPaX<4) 

PHIPaX(2) 

CALL  TEMP  IPSItTHTAiPHItPSlP*THTAPtPHIP*O.D«OtO. 0*0 t0.D«O«O.O«Ot 
IO.O*OtO*l)*OfK) 

call  MBAR|MXOIXtMYOI»MZOItPSI*THTAtPHIfPSlP*THTAPtPHIP) 

DERXCll  a PHIP 

DERX<2)  » MXOIX  ♦ PSIP*THTAP«OCOS(THTA)*MZOI*DTAN(THTA) 
1*(2,D*0-HATIO)»PS1P*THTAP»DTAN<THTA)*OSINITHTA> 
2*RATI0*PMIP*THTAP*0T AN { THTA) 

DERXOl  a TMTAP 

DERX(AI  a MYOI  -« 1 .0*0-RATI0» • «PSIP»*2) * (OSIN (THTA » AOCOS (THTA> ) 
I-RaTI0*PHIP*PS1P*0C0S(THTA) 

DERX(5)  » PSIP 

DERX(6>  a MZOI/OCOS(THTA) • (2.D«0>RATIO)*PSIP*THTaP* 

IDTAN(THTa)  • RAT10»PHIP*THTAP/'0C0S(THTA» 
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NEQ»(NC*U*6 
DO  5 !B7tNEQ 
S DERX(I>sO.O«0 
DO  10  l»ltNC 

1F(NC0N(1) .E0*0>  GO  TO  10 
DERX(6*1)  » X(6«NC«n 
DEHX{2*NC*6*I>  ■ X(3*NC*6*I) 

0ERX(4*NC«6*I>  « X<5*NC*6*1) 

PI  > X(6*n 
P2  B X(2*NC*6«1) 

P3  B X(4*NC«64l) 

P1P»  X(NC*6«1) 

P2P»  X(3»NC*6*I) 

P3PB  X(5*NC*6*I) 

call  TEMP<PSltTHTA.PHIfPSIP.THTAP.PHIP.Pl.P2.P3.PlP»P2P.P3PtK) 

CAUL  MBARA(Ml*M2fM3»PSI»TMTAiPHItPSlP»THTAPiPHlPtPl#P2»P3»PlP»P2Pt 
lP3Ptl) 

OERX (6*NC* I ) bMI ♦THTAP*P3P*OCOS (THTA) ♦PSIP*P2P*DC0S (THTA) 
1-PS1P*THTAP*P2*0SIN(THTA>  ♦ M3*DTAN (THTA)  ♦ 
2MZ0l*P2/(0C0S(THTA>B*2)  ♦ (2.D*0-RATIO)*OTAN(THTA>* 
3DSIN(THTA)b(THTAP»P3P  * PSIP*P2PI ♦ (2«0*0-RATIO) ♦PSIP^THTAP 
4*P2«DSIN(THTA)*(1.D*0*1.0*0/(DCOS(THTA)**2)  ) * 
5RATI0*(P1P*THTAP*DTAM<THTA)>PHIPBP2P»0TAN(THTA) 

6*PhIP*THTAPbP2/ (0C0S<  THTA) **2) ) 

OERX (6*3*NC*I ) ■M2» ( 1 ,D*0»RAT10) •2.D*0*PSIP*P3P*OSIN (THTA) 
1*DCOS(THTA)-(1.O*0-RATIO)bP2»(I.D*O-2#O*0*(OS1NITHTA)**2) )* 
2PSIP«*2  - RAT10B(P1P*PSIP*DC0S(THTA)*PH1P»P3P*DC0S(THTA) 
3.PHIP»PS1P«P2*DS1N(THTA) ) 

DERX (6*S»NC*I ) BMS/DCOS ( THTA) *MZ0I*P2*0TAN (THTA) /OCOS (THTA) 
14(2.D*0-RATIO)*OTAN(THTA)*(PS1P*P2P«THTAP»P3P)  ♦ 

2 (2.0*0-RAT 10) •PS1P*THTAP»P2/DC0S (THTA) BBZ  ♦ 

3RAT10*P1P*THTAP/DC0S (THTA) ♦RAT  10* (PH1P*THTAP*P2B0TAN( THTA) 

A/OCOS ( THT A ) ♦PH IP*P2P/DC0S ( THT A ) ) 

10  continue 

RETURN 

END 


SUBROUTINE  PART(SAMPLEtNOELTA»NtNTOT«NSTART) 

C 

IMPLICIT  REAL*8  (A-H»0-Z) 

COMMON/MEA5/PHIM(1000)  tTHTAM(lOOO)  (PSIMdOaO) 
COMMON/SOL/T(1000) iPHI (1000) tTHTA( 1000) tPSI (1000) 
00  13  la2tNT0T 
IF(T(1) .NE.0.0^0)  NTEHPbI 
13  continue 

TMAXbT(NTEMP) 

00  11  laNSTARTtNTOT 
Jb(I-1)*N0ELTA^1 
1F(J.6T*NTEHP)  GO  TO  11 

RATIOb(T(J)-T(NSTART) )/TMAX 
1F(RAT10*6T. SAMPLE)  GO  TO  11 
NbI 

T(1)bT(J) 

PH1M(I)bPHIM(J) 
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THTAM(I»«THTAM(J) 

PSIM(I)aPSlM(J) 

11  continue 
NaN*l-NSTART 
DO  12  I"ltN 
JaMNSTAHT-l 
T<I)»T(J) 

PHIM(1)»PHIM(J> 

THTAM(I)aTHTAMIJ) 

12  PS1H(1)>PS1H(J) 

NlsN*l 

IFtNTOT.LT.Nl)  60  TO  14 
DO  15  laNltNTOT 
15  T<1)>0.0*0 

14  continue 

MRlTE(06«10)SAHPLEfNDELTA*NSTART*N*T(l) «T(N) 

10  FORMAT(lMltlX*»»*»»  DATA  BLOCK  ALTERED  •••*• */t 1 Xi • SAMPLE 
lF7.4,/,lAt*N0ELTA  * I5*/t 1 X« *NST ART  >• i ISt/t IXi *N  ■**I5t/t 

21X**T(1)  ■ • fF7.4,5X»'T (N)  ■ ♦tF7,4./# 

31Xt2B('«'l > 

RETURN 

END 


FUNCTION  OUM(ItJ) 

IMPLICIT  REAL*e(A«HtO-Z) 
DUMa0«D*0 

IFd.EOtJ)  DUMal.O«0 

RETURN 

END 


SUBROUTINE  ICSIACON) 

IMPLICIT  REALMS  (A-H«0-Zl 
DIMENSION  ACONOO) 

COMMON/MEAS/PHIM(1000)  (THTAMdOOO)  *PSIM(1000) 
COMMON/SOL/T IIOOO) fPHi 11000) tTHTAl 1000) tPSl (1000) 
DTaT(2)-T(l) 

PHlOaPHlM(l) 

THTAOaTMTAM(l) 

PSlOaPSlM(l) 

AC0N(25)aPHI0 

ACON(27)aTHTAO 

AC0N(29)»PS10 

AC0Nl26)a(PHIM(2>-PHI0)/DT 

ACON (28) a (THTAM(2) -THTAO) /DT 

ACON(30)a(PSlM(2)«PSIO)/UT 

RETURN 

END 


C 


SUBROUTINE  INV(C*NCtNCSltNXX*CC) 
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IMPLICIT  ReAL*8(A-H,0-Z) 

DIMENSION  C(30*31)  tMXXOOl  tPlV0T<2l  fCC(30i60) 

3000  FORMAT</IOX««DET  is  equal  to  ZERO') 

NCT*NC*2 
NCPl>NC«I 
00  10  I«ltNC 
00  10  jaltNC 
10  CC<lt J)»C(ItJ> 

00  20  laltNC 
00  20  JaNCPlfNCT 
20  CC(I«J)«0.0«0 
00  30  1>1«NC 
30  CCIItNC*l)«l.D«0 
DO  205  l>ltNC 
PlVOT(l)»CC(Itl) 

00  200  KaltNC 
PIV0T(2)*CC(K*I) 

126  IF(K-I)  13S*130«140 
130  DO  150  JaltNCT 

1F(P1V0T(1))  13Af210tl34 

134  CC(Kf J>«CC(ItJI/PlVOT(l) 

ISO  continue 

GO  TO  200 

135  00  160  UaltNCT 
1F(P1V0T(1)>  136*160tl36 

136  CC(K,J)«CC(K«J)-CC(If J)*P1V0T(2)/PIV0T(1) 

160  continue 

GO  TO  200 
140  00  170  JaliNCT 

IF(PIV0T(2))  14Stl70fl4S 
145  CC(K,J)«CC(K«J) /PIVOT {2>-CC(I«J) 

170  continue 
200  CONTINUE 
205  continue 
GO  TO  250 
210  WRITE  (6*3000) 

250  00  300  laltNC 
DO  300  JaltNC 
300  C(ltU)aCC(I*U*NC) 

NCS«NCSl 

350  lF(NCS-NC)500«500t400 
400  00  420  laltNC 
WXX(I)aC(l»NCS) 

420  C(ItNCS)«0.0*0 
DO  450  I"1«NC 
00  450  JaltNC 

450  C(lfNCS)«C(I«NCS)«C(ItJ)*WXX(J) 

NCS»NCS-1 
GO  TO  350 
500  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  NOISE (SPHl tSTHTAfSPSI tANPHl « AMTHTA* AMPSI tN* 
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ISAHPLEtNDELTAtNTOT) 

THIS  SUBHOUTINE  ADOS  GAUSSIAN  NOISE  TO  DATA 
IMPLICIT  REAL*8U-H*0-Z> 

COMMON/MEAS/PHIMIIOOO) (THTAMIIOOO) tPSIM(lOOO) 
COMMON/SOL/T(lO00l  tPHKlOOO)  tTHTA(lOOO)  tPSl(lOOO) 

MRlTE(06f 100)SPSliSTHTAtSPHI,AMPSI*AHTHTA«AMPHlfN 
call  TIMEIlYEARtlSEC) 

7 IF(ISEC-IOOO)  5«5t6 
6 ISEC>ISEC-1000 
60  TO  7 

S lXa<lSEC«IYEAR)*2*l 
IYsIX*2 
IZ*1X«4 
SlGlPHsOtO^O 
SIG2PHa0.0«0 
SI63PHxO.O*0 
SIGlTHaO.O«0 
SlG2THaO.O«0 
SIG3THa0.0«0 
SIGlPSa0.0«0 
SI62PS«0.D*0 
SI63PSa0.D«0 
TOTPHIaO.D^O 
TOTTHAaO.D»0 
TOTPSIaO.D*0 
LOOPal 
NLOOPal 
00  lO  lal.N 

call  GAUSSdXtSPHItAMPHltVPHIl 
call  gauss (lYiSTHTA,AMTMTAfVTHTA> 

CALL  GAUSS(lZfSPSI»AHPSl«VPSI) 

PHlM<I)aPHlH(I)*VPHI 

THTAM( I ) alHTAM ( I ) * VTHTA 

PSIM<naPSIM(I)»VPSI 

TOTPHlaTOTPHl*VPHI 

TOTTMAaTOTTMA*VTHTA 

TOTPSIaTOTPSl*YPSI 

PaOABS(VPHl) 

QaOABS(VTHTA) 

HaDABS(VPSl) 

IF(P.LT.SPMl)SIGlPHaSI6lPM*l,D*0 

lF(P.6E«SPHI.ANO*P*LT.2tO«0*SPHl)SIG2PHaSlG2PH«l.D«0 

lF(P.6E.2,0*0*SPHI«ANO.P.LE.3.D«0*SPHI>SI63PHaSIG3PH«1.0«0 

lF(O.LT.STHTA>SI61THaSI6ITH«l.D«0 

IF(Q.6E.STHTA.ANO.O.LT.2.O«0*STHTA)SIG2THaS162TH»l.O*0 

lF(Q.6E>2.D«0«STHTA.AND.Q.LE.3tD«0*STHTA)SlG3THaSIG3TH*l.D*0 

IF(R.LT.SPSI)SIGlPSaSIGlPS*l.D*0 

lF<R.6E.SPSI.AND>R.LT.2.0*0*SPSI)SIG2PSxSI62PS*1.0«0 

lF(R.6E.2.0»0*SPSl*AND.R.LE.3.0«0*SPSl)SIG3PSaS163PS«l.D*0 

10  continue 

WRITE (06* 103) TOTPSItTOTTHAtTOTPHltSIGlPStSIGlTHiSIGlPHf 
lSIG2PS*S162THtSIG2PHfSlG3PStSIGaTHtSI63PH 
100  FORMATIlHltlXt****  GAUSSIAN  NOISE  ADDED  •**'«/tlXt 

I'ONE  standard  DEVIATION  - PSI  ai ,F10.6*5Xf 'THETA  ai,Fl0.6i5Xt 
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2'PHl  ••*F10.6*/tlXt'CONSTANT  BIAS  • PSI  ■•tFlOtBtSX* 

3'THETA  ■•tF10«6t5Xt'PHl  «• tFl0.6t/*lXt 'TOTAL  NO  OF  TINE  POINTS  USE 
AO  IN  FIT  (N)a*»IS) 

103  FORMATdNlflX.'SUM  OF  NOISE  ADDED' #/#6X* 'PSI  • i6Xt 'THTA' •TX* 'PHI » , 
I/tlX«3Flt).3i//*lX«'N0ISE  CONTRIBUTION  WITHIN  ONE  SIGHA't/t 
26X • 'PSI * *6X •' THTA ' t7X •' PHI '•/»lX»3F10.3t//* 

31Xt 'NOISE  CONTRIBUTION  BETWEEN  ONE  SI6MA  AND  TWO  S16HA't/t 
A6Xf 'PSI' «6X(' THTA' fTXt 'PHI' •/* IX t3F10»3t//« IX t 'NOISE  CONTRIBUTION 
5BETWEEN  TWO  SIGMA  AND  THREE  SIGMA ' »/t6Xt 'PSI ' t6X» 'THTA' tTX# 'PHI ' t 
6/.1X.3F10.3) 

RETURN 

END 

C 

C 

subroutine  RANDUdXiIYtYPL) 

C 

IMPLICIT  REAL*8(A-HfO>Z> 

IY>IX*65S39 

IF(lY)5t6«6 

5 lV»lY*21A7Ae3647*l 

6 YPL“IV 

YPL"YPL»*A656613D-9 

RETURN 

END 


subroutine  6AUSS(IX«S*AM»V) 

IMPLICIT  REAL*6(A*HfO-Z) 

AaOiD«0 

DO  SO  Ial*A8 

call  RANDU(IXtIY»Y) 

IXalY 
50  A«A«Y 

Vs(A-2A«D*0)/2«0*0 

VaV*S*AM 

RETURN 

END 


SUBROUTINE  INT6RL<FCT*DRK6S(NtERR0R*AUXtYt0ERYt0TMIN) 

IMPLICIT  REAL«8(A-Ht0*Z) 

DIMENSION  AUX(6«1I tY(l) iDERY(l) 

COMMON/MEAS/PHlMdOOO)  *THTAM(I000)  tPSIMdOOO) 
COMMON/SOL/T<  10001  tPHI  (1000)  tTHTAd 000)  fPSI  dOOO) 
COMMON/COEF/CON(30) tCCON(30)tNC«NCON(30) 
common/m A TRIX/C 130 160 ) t RESSUM t SPHI • STHA  t SPSI 
A>6.28318S310«0 
NL>NC-S 
DO  10  l>lt6 
10  Y(I)aCONINL*I-l) 

DO  11  laltNC 
Y(6*naDUM(NL*l) 

DO  11  J«ltS 
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11  Y(6*J*NC«l)B0UM(NL*JtI> 

NEQ>(NC«1>*6 

DNEQaNEQ 

NL>NC*1 

K»1 

PHI<K)sY(l) 

THTA(K)»Y(3) 

PS1(K|«Y(5> 

THTA(K>»THTA(K>*C0N(23) 

PSI(KI>PS1(K)«C0N(24) 

DPHlaPHlM(K)-PHI (K) 

1F(UPHI.6T.<A/2«D*0))  PHI (K) >PHI (K) *A 
IF(OPH1,LT.I-A/2.0*0) I PHI (K)-PHI(K)-A 
00  19  I»lfNC 

Y(6*2*NC«I)«Y(6«2«NC*II *0UM(23*I) 

Y<6*4«NC«I)aYI6«4*NC*I)«0UM(24tn 

19  continue 

00  12  lal*NC 
00  12  J>ltNC 

C(I*JI«Y(6*I)*Y(6*J)/SPH1»*2  ♦ Y<6*2*NC*I)*Y(6*2*NC*J)/STHA**2 
!♦  Y<b*4*NC*l)»Y(6*4*NC*JI/SPSI**2 

12  C(1*J)«C(1»J)*CC0N(I)»CC0N(J) 

00  17  J>ltNC 

C<J»NL»«(PH1M<K>-PHI(K) »*Y(6*J)/SPHI*»2  ♦ < THTAM (K) -THTA (K>  » * 
1Y<6*2»NC*J)/STHA»*2  ♦ (PSIM(K>-PSI IK) >*Y(6*4*NC*J)/SPS1**2 
17  C(JtNLI"C(JtNLI*CCON(J) 

RESSUM  » ((PHIM|K)-PHI(K)>/SPHI)*#2  ♦ I ( THTAM(K) -TMTA (K) ) /STHA) **2 
1*  (<PSIM(K)-PS1(K))/SPSI)»*2 
00  13  Ks2tN 
XIaT (K»l) 

0TE*T(K)-T(K-l> 

DTaOTE 
NOT"  I 

1FI0TMIN.6E.0TE)  60  TO  18 

not»ote/utmin 

OTaNOT 

OTaOTE/OT 

16  continue 

call  ORKGSIXltYfOERY*NEQ«FCT«AUX»DT«NDT) 

IFIYIl) .LE. (•A>)Y(l)aY(l)*A 
IF(Y(l>.GE*A)Y(l)aY(I)-A 
PHI iK)aY(l) 

THTA(K)aY(3) 

PSl  lK)aY(5) 

THTA(K)aTHTA(K)*C0N(23) 

PSl (K)8PSI (K}*C0N(24) 

0PHl8PHlM(K»-PHI (K) 

IFIOPHI.LT* (•A/2. 0*0))  PHI(K)aPHl (K)-A 
IF (OPHI. GT. (A/2.0.0) ) PHl(K>aPHI (K)*A 
DO  20  lal.NC 

Y(6*2»NC.I)aY(6*2*NC«I)*0UM(23,I) 

YI6*4*NC*l)aY(6.4*NC.I) *0UM(24.I) 

20  continue 

DO  IS  lal.NC 
00  IS  jal.NC 
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C(ItJ)«Yt6*I)*YI6*J)/SPHl**2  • Y(6*2»NC*I)«Y<6*2*NC*J> /STHA«»2 
!♦  y(6*4«NC*II»YI6^4*NC*JI/SPS1««2  ♦ C(IiJ) 

15  C(ltJ)sC(ltJ)*CCON(I)*CCON(J) 

00  15  JalfNC 

C(JtNL)»IPHlM(K)-PHI(K) )*Y(6*J)/SPHI»*2  ♦ (THTAM(K)-THTA(K) )♦ 
1Y(6*2*NC*J>/STHA**2  ♦ <PSIM(K)-PS1 (K) ) *Y (6*4*NC>J> /SPS1**2*C ( JtNLl 

16  C(JtNL)«C(JfNL)*CCON(J) 

RESSUM  a ( (PHIM IK) -PHI (K) ) /SPHI ) **2  » I (THTAM |K> -THTA (K) ) /STHA) aaZ 
I*  ((PS1M(K)-PS1(K))/SPSI>**2  ♦ RESSUM 
13  continue 
RETURN 
END 


SUBROUTINE  TEMP (PSI «THTA*PHI tPSIPtTHTAPtPHIPi PI (PBtPatPlPtPBPtPaPi 
IK) 

IMPLICIT  REAL*8(A-HfO-ZI 
REAL*8  MAG 

C0MM0N/TC0EFl/SPN«AtB«MA6tQA*RAtAD*B0«ALPAiBTA»ALPAP«BTAP«DLTA2*Pf 

iQtRfSMAGtOLTA 

C0MM0N/C0EF/C0N(30) *CCON(30) «NCiNCON<30l 
SPN"PIP-PSIP*P2*0C0S ITHTA) -P3P*OSIN (THTA) 

A«-P3*DS1N (PSI I *05IN(THTA) «P2*OCOS (PSI )*OCOS ITHTA) 

Ba-P3*DC0S(PSI) 

MAG*  2,0>0*(P2aDSIN(THTAI*DCOS(PSI> 

1*  P3*  0C0S(THTA)*0SIN(PS1) )«DCOS(PSI)aDCOS<THTA) 

QA  ■ P2P 

RA  ■ P3P«DC0S<TMTA)-PSIP*P2*DSIN<THTA) 

AO—P3P«OSIN  (PSI ) aOSlN  (THTA)-P2*  <PSIP*DS1N  (PSI ) •OCOS  (THTA) 

1 ♦THTAPaOCOS (PSI ) aOSIN (THTA) ) -P3» (PSIP*DCOS (PSI ) •OSIN (THTA) 
2«THTAP*0SIN (PSI ) *OCOS (THTA) ) «P2P«  OCOS (PSD aOCOS (THTA) 

B0»-  P3P*DC0S ( PS 1 ) ♦ PS I P*P3*DS I N ( PS I ) 

ALPA  > OCOSIPSI)  • DSIN(THTA) 

BTA  >-DSII^(PSl) 

ALPAP  > -PSIP*0S1N(PS1)*0SIN(THTA)«THTAP*0C0S(PSI)*DC0S(THTA) 

BTAP  • -PSIP*0C0S(PS1) 

0LTA2  • (OCOS(PSI)*DSINITHTA) )«*2  * (DSIN(PSI) ) **2 
P ■ PHIP  - PSIP  • OSIN(THTA) 

Q > THTAP 

R ■ PSIP*OCOS(THTA) 

if(k.eq*0)  go  to  10 

ALPAaTHTA 

BTAaOATAN (-DTAN (PSI ) /OCOS (THTA) ) 

Aap2 

Ba. (P3/ (OCOS (PSD ••2) «DTAN (PSI ) *OTAN (THTA) aPZ) • 
l(DC0S(BTA)*a2)/0C0S(THTA) 

0LTA2aALPAaa2»BTA**2 
MAGa2.0*0* (ALPA«A«BTA*B) 

10  continue 

0LTA«0SQRT(0LTA2) 

SMAGbO.0«0 

IF(0LTA2.EQ. 0.0*0)  GO  TO  9 
SMAGaMA6/(2.D*0*OLTA) 

9 continue 
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A0b0.D*0 

B0>0,D«0 

ALPAPaOtO^O 

eTAPaO.0^0 

RETURN 

END 


SUBROUTINE  RE60UT (ACONtWXX) 

IMPLICIT  REAL*a(A-H*0-Z) 

REAUaS  MPRMtIPRM 

OIMENSION  ACONI30)  tWXXOO) 

COMMON/COEF/CON (30) t CCON (30) tNCtNCON (30) 

COMHON/C0NST/RAT10«MPRM«lPHH«Af tangle 

Cl  a IPRH/MPRM 

DO  10  I«it22 

ACON(I)  a MPRH«C0N(I) 

NXXd)  a MPRMaWXX(l) 

ACON(l)  a AC0N(1)*C1 

AC0N(3)aAC0N(3)*Cl 

AC0N(5)  a AC0N(5)*C1 

MXXd)  a MXXd)  • Cl 

wXX(3)aWXX(3)aCl 

MXXlS)  a nXX(5)  aCl 

AC0N(21)aAC0N(21)aCl 

AC0N(22)aAC0N(22)aCl 

MXX(21)aMXX(21)*Cl 

WXX(22)a»IXX(22)*Cl 

ACON(10)aACONdO)aCiaRATlO 

ACON (11) aACON (11) aRAT I 0 

AC0N(13)  a AC0N(13)  a RATIO 

AC0N(14>  a ACON(U)  a RATIO 

WXXdO)a«IXXdO)aCiaRAT10 

WXX(ll)a«IXX<ll)aRAT10 

WXX(13)  a WXX(13)  a RATIO 

WXX(U)  a inxxdA)  a RATIO 

DO  11  Ia23,30 

WXX(I)aiKAX(I)/A 

ACON(I)  a CON(l)/A 

RETURN 

END 


SUBROUTINE  REG (ACON) 

IMPLICIT  REAL*8(A-H«0-Z) 

REALas  MPRM*1PRM 
DIMENSION  ACON(30) 

COMHON/COEF/CON ( 30 ) vCCON ( 30 ) « NC  *NCON ( 30 ) 

COMMON/CONST/RATlOiHPRMdPRMt  At  TANGLE 

DO  lO  lal,22 

CONd)  a ACON(l)/MPRM 

CONd)  a ACON(l)/IPRM 

C0N(3)aAC0N(3)/lPRM 
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C0N(5)  ■ AC0N<5)/IPRH 
CON ( 1 0 ) aACON  < 1 0 ) / ( IPRM*RAT 10) 
CON  (ll)aCON(in /RATIO 
CONUS)  a C0NI13)/RATI0 
C0NU4)  a CONUA)/RATIO 
C0N<21) aACON (21 )/IPRM 
CON  < 22 ) a ACON ( 22 ) / 1 PRM 
00  11  Ia23*30 
II  CONU)  a AC0NII)*A 
RETURN 
END  . 


SUBROUTINE  MBARA(MltM2»M3tPSI*THTA*PHI»PSlPtTHTAP»PHlP* 
lPltP2«P3tPlPtP2PtP3Pf J) 

IMPLICIT  REAL*S(A-HtO-Z> 

REAL«8  Ml*M2tM3»MA6tMYTtMZT«MYTPtMZTPtHYTE«MZTE 
COMMON/TCOEFl/SPNtAtBfHAG«QAtRAtAO*BO*ALPAtBTA» 
lALPAP.8TAP(0LTA2«PtQtR*SMA6tDLTA 
COMMON/COEF/CON(30)fCCON(30) »NCtNCON(30> 
C0HM0N/C0NST/RATI0tZliZ2tZ3tI4 
MlaSPN*(CON(13)«CONUA)aOLTA2l*P*(OUMI13tJ) 
l*DUM(14tJ)aOLTA2  * C0N(14)*MAG) 

HlaHl ♦DUM ( 1 0 1 J>  «SPNaCON (11) aOLTA^PaOUH ( 1 1 1 J) *DLTA*P*CON (11) aSHAG 
TEMPI  a CON(l)  « C0NI5)«0LTA2 
TEMP2  a UUM(l»J>«0UM(StJ)a0LTA2  * CQN(5)*MAG 
TEMP3Y  a CON (2) *CON (6) »0LTA2-C0N (7) aBTAaaZ 
TEMP3Z  a C0N(2)  ♦ C0NI6)»DLTA2  - C0N(7) *ALPA**2 
TEMP4Y  a DUM(2«J)«^DUM(6iJ)*0LTA2  ♦ C0N(6>*MA6 
1-OUM ( 7 1 J) aBTAaaZ-Z.OaOaBTAaCON ( 7) aB 
TEMP42  a 0UM(2*J)  » OUM (6* J) aDLTA2  * C0N(6)aMA6 
l-0UH(7tJ)aALPAa«2  . 2.D*0aALPAaC0N (7) aA 
TEMPlaTeMPl*C0N(3)aDLTA 
TEMP2aTEMP2«DUM(3f J)aDLTA»C0N(3)aSMAG 
TEMP3YaTEMP3Y*CON(9)aOLTA 
TEMP3ZaTEMP3ZaC0N (9) aQLTA 
TEMP4YaTEMP4Y«0UM (9tU) a0LTA«C0N (9) aSMAG 
TEMP4ZaTEMP4ZaOUM (9« J) aOLTA4CON (9) aSMAG 
TEMP5Y  a CON(3)*CON(9)aOLTA2-CON(10)aBTAa*2 

TEMP5Z  a C0N(3)  ♦ C0N(9)a0LTA2  - CON(10)aALPAa#2 
TCMP6Y  a DUMOf J)*0UM(9tJ)a0LTA2  « C0N(9)aMAG 
l-OUM(10«J)aBTAa*2  . 2.0«0aC0N 1 10) aBTAaB 
TEMP6Z  a OUM ( 3t J) *OUM (9* J) aOLTAZ  * C0N(9)aMA6 

l-0UM(10iJ)aALPAaa2  . 2.D40acON t 10) aALPAaA 

TEMP7  a C0NI4)*C0N(12)a0LTA2 

TEMPS  a UUM<4*J)aOUM(12tO)aDLTA2  * C0N(12)aMAG 
M2  a AaTEMPl  ♦ ALPAaTEMP2  ♦ QAaTEMPSY  ♦ QaTEMP4Y 
l40UM(8tJ)aRaALPAaBTA4C0N(8)a(RAaALPAa6TA4RaAaBTA 
2«RaALPAaB) 4AOaTEMPSY  ♦ ALPAPaTEMP6Y 

340UM(ll«J)aBTAPaALPAaRTA  a CON(ll)a(BOaALPAaBTA«BTAPaAaBTA 
4*BTApaALPAaB)  ♦ SPN*BTA  ajEMPT  ♦ BapajEMPT 
54paBTA  atEMPS 

M3a-8aTEMPl-BTAaTEMP24RAaTEMP3Z  ♦ RaTEMP4Z 

1 *0UM (8t0) aQaALPAaBTAaCON (8) a lOAaALPAaBTAaQaAaBTA 
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c 

c 

c 


2*Q*ALPA*a)-B0*TEMP5Z  - 8TAP*TEMP6Z 

3-DUM ( 1 1 1 J) •ALPAP*ALPA*BTA-CON  < 1 1 ) • ( AD*ALPA*BTA»ALPAP*A*BTA 
♦ ♦ALPAP*Al.PA«B>*SPN*ALPA*TEMP7  ♦ A*P*TEMP7 
5*P«ALPA*TEMP8 
MYT»C0N(2n*THTA 
MYT»MYT*C0N (20) *THTAP 
MZT«C0N(22)*PSI 
MZT»MZT*CON (19) •PSIP 
MYTP>OUM (21 t J) •THTA*CON (21 ) *P2 
MYTP«MYTP«0UM (20f J) •THTAP*C0N (20) •P2P 
MZTP»0UM(22tJ)*PSl*C0N(22)»P3 
MZTP»MZTP*0UM(19»J)*PSIP*C0N(19)»P3P 
MYTE» (OUM (21 1 J) *THTA*CON (21 ) •PZ) /RATIO 
MYTEBMYrE*(0UH(20f J)*THTAP*CON (20) OPZP) /RATIO 
MZTEb (DUM (22* J) •PSI •con (22) •PB) /RATIO 
HZTEbMZTE«(DUM(19(J)«PSIP*C0N(19)»P3P)/RATI0 
M1«M1 ♦MYTE*DSIN (PSI ) •OCOS (THTA) ♦MYT» (P3«DCOS (PSD  *DCOS ITHTA) 

1 -P2*0SIN(PSI)*DSIN(THTA) )/RATI0 

2 -MZTE*OSIN ( THT A ) -MZT*P2*0C0S (THTA ) /RATIO 
M2BM2«HYrP»0C0S (PSI ) -MYT«P3«DSIN (PSI ) 

M3-M3*mv  TPOQSIN (PSI ) •OSIN ( THTA) ♦MYT^ (PS^OCOS (PSI ) •OSIN ( THTA) 

1 ♦PZ^OSIN (PSI )^DCOS (THTA) ) 

2 •MZTP^OCOS (THTA I -MZT^PZ^OSIN (THTA) 

RETURN 

ENO 


SUBROUTINE  MBAR(MXOlXfMYOI*MZOI»PSItTHTAtPHI«PSIPtTHTAPfPHIP) 

IMPLICIT  REAL^8(A-HtO-Z) 

REALMS  MX0IXtMY01tMZ0I»MA6fMYTtHZT 

C0MM0N/TC0EFl/SPNtA*BtMA6t0AiRAtA0«BD«ALPA«BTA*ALPAPtBTAP<0LTA2t 

lPtQ*R«SMAGtOLTA 

COMMON/COEF/CON(30I  tCC0N(30|  iNCtNCONOO) 

C0MM0N/C0NST/RATI0»Zl*Z2tZ3tIA 

HX0lX>  P^(C0N(13)*C0N(1A)^0LTA2) 

MYOI  a ALPA^(CON(1)«CON(S)^DLTA2)«Q^(CON(2)«CON(6)*DLTA2 
l-CON(7)^ijTA**2)  •C0N(8)*R^ALPA*BTA*ALPAP*(C0N(3)  ♦ 

2C0N (9) •0LTA2-C0N ( 10 ) •BTA^^Zl •CON (11) •BTAP^ALPA*BT A 
3^P^BTA  •(C0N(4)^C0N(12)^0LTA2( 

MZOI  ■ -dTA^(CON(ll •C0N(5)^0LTA2)^R^(C0N(2)^C0N(6)«0LTA2 
1 -CON ( 7) •ALPA**2) •CON (8) •U^ALPA*BTA-BTAP^ (CON (3) • 

2CON(9)«OLTA2  - CON( 1 0 ) •ALPA**2) -CON ( 1 1 ) •ALPAP^ALPA*BTA 
3^P*ALPA^ (CON (4) •CON( 12) •OLTAZ) 

MXOlXaHXOlX^CON ( 10) •P^CON (11) ^OLTA 

MYOlaMYODALPA^CONO)  •0LTA^Q^C0N(9)  •OLTA 

MZ0l3MZ0I-BTA^C0N(3)^0LTA^R^C0N(9)^0LTA 

MYT»C0N(21I^THTA 

MYT«MYT^CON(20>^THTAP 

MZT»C0N(22»«PSI 

MZI»MZT^C0N(19)^PS1P 

MXOlXaMXOIX^(MYT^DSlN(PSl)^OCOS(THTA)-MZT^OSIN(THTA))/RATIO 

MYOIbMYOI^MYT^OCOS(PSI) 

MZ01«MZ01^MYT^0S1N{PSD^0SIN(THTA)^MZT^0C0S(THTA) 

retuhn 

ENO 
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Table  2.  Fortran  Listing:  Asymmetric  Parameter  Extraction  Program 


A N 6 L E S 

VARNER  VERSION  1 APR.  1977 

COMPUTES  AERO  COEFFICIENTS  FOR  MIRROR  SYMMETRY  BODY 
BODY  AXIS  SYSTEM.  AERODYNAMIC  ASYMMETRY 
EXACT  3 jOF  EQUATIONS  UTILIZED 


IMPLICIT  REAL*8(A-H.0-Z» 

REAL*8  IX.IY.IZ.IZXf JX.JY.JZ.IXP.IYP.IZP. MXlfMX2iMYl.MY2.MZl.MZ2 
REAL*e  AL1ST{72) 

REAL*8  CARD  1 10) 

REAL«A  TS(IOOO) .PHIS (1000) . THTASd 000 ) tPSIS (1000) 

DIMENSION  AUX ( 5.438) ♦Y(438) . DERY (438) .WXX (72) .ACON(72) tCT (72.73) 
l.WX(72) .AAC0N(72) fCC(72.l44) 
external  FCT.DRKGS 

COMM0N/MEAS/PHIM(1000)  .THTAMdOOO)  .PSIMdOOO) 

COMMON/SOL/TdOOO)  fPHl  (1000)  .THTAdOOO)  fPSI  (1000) 
COMMON/COEF/CON(72) .CCON(72) .NC.NCON(72) 

COMMON/MATRIX/C (72.73) *RESSUM»S16PHl.SIGTHAt  S16PS1 
COMMON/CONST/ JX  fJY.JZ. iXP. I YP. IZP. I ANGLE . IX. lY. IZ 
COMMON/CONST2/MXlfMX2.HYlfMY2fMZlfMZ2»A 

equivalence  (AUX(1)  tCTd) ) . (Yd)  »CT(2191) ) . (OERY ( 1 ) tCT (2629) ) t 
l(CCd).C(l))  ^ 

• (PHISd)  fPHId) ) « (THTASd)  .THTA(I)  ) . (PSISd ) tPSI  ( 1 ) ) 

. 


DATA  ALIST/'CM  0 
1*CM  A3  *.»CM  A2 
2*CH  MAB  '.'CM  A2 


• f 'CM 
MA'.'CM 
MB'. 'CM 
'.'CMQ 


A 

MB 

MAB 


3'CMQ  A '.'CMQ  MA  '.'CMQ  A2 

4'CMQ  A MU'. 'CMQ  MAB  '.'CMR  B 

5'CN  A B '.'CN  MA  B '.'CN  B3 
6'CNR  0 '.'CNR  A '.'CNR  MA 

7'CNR  A MAt.'CNR  A MB'. 'CNR  MAB 


8*CL  A 
9'CLP  MA 
••IZX 

••floang 

«'R(0) 


'.'CL  B 
'.'CLP  MB 
'.'CNR  TUN 
P'.'PHKO) 
'/ 


'.'CL  A B 
'.'CLP  A2 
'.'CMQ  TUN 
'.'THTA)0) 


.'CM  MA 

'.'CM  A2 

'.'CM  A MA 

.'CM  B2 

'.'CM  82  MB 

'•'CM  A MB 

.'CM  A B2 

'.'CM  MA  B2 

'.'CMQ  0 

f'CMQ  A MA'.'CMQ  MB 

'.'CMQ  B2 

.'CMP  B 

•f'CN  0 

'.'CN  B 

.'CN  A2  B 

'.'CN  B MB 

'.'CN  B MAB'. 

.'CNR  A2 

'.'CNR  MB 

'.'CNR  B2 

f 

.'CNQ  B 

'•'CNP  A 

'.'CL  0 

1 f 

.'CL  MA  B 

'.'CL  B MB 

'.'CL  P 

• t 

•'CLP  82 

*. 'DUMMY 

'.'DUMMY 

f f 

.'Cma  tun 

'.'CNB  TUN 

'.'FLOANG 

T*# 

.•PSKO) 

'•'P(0) 

'.•Q(0) 

1 1 

NAMELIST/NAM/NSWTCH.NTOT.NC.NITER.OTMIN.HOUND.AREA.OIA.IX.IY.RO.U. 
lAACON.NCON.NSESWH.SPHI.STHTA.SPSl.AMPHI.AMTHTA.AMPSI. 
2SAMPLE.N0ELTA. lUPDAT. IZ. IZX . I angle. INLCON. OTGEN 
3. IGPNO. IBLKNO. IRCGEN.SIGPHI .SI6THA .SIGPSI .NSTART. IRCSTR 
100  format (1H1.20X. '*•**••*•  INPUT  DATA 

UX. 'GROUP  NUMBER  ■'. 15. 5X. ' BLOCK  NUMBER  »'.I5./. 
19X.'TIME'.10X»'PSI'.10X.'THETA'.10X.'PHI'.//) 


101  FORMAT(7F10.0) 

102  F0RMAT(5<.4(E13.6.1X)). 

103  F0RMAT(5X.2(A8.2X.2(E13.6.5X)) ) 

104  F0RMAT(/. IX. 'ITERATION  NO' . I5./.2 (21X. 'CON' . 1 IX. 'DELTA  CON')) 

105  F0RMAT(2(11X. 'FINAL  CONSTANTS  AND  PROBABLE  ERROR' . lOX) ./) 

106  F0RMAT(2(5X.iC('.13.')»'.E13,6.SX.' ('.E13.6.') '.2X.A8)) 

107  FORMATdHl.SX.'DATA  FIT  COMPARISON' .//.IIX.  'TIME'. 
13X.'PSI/FIT'.3X,'PSI/EXP'.3X.'DEL  PSI ' .2X. 'THT A/FIT ' .2X. 'THTA/EXP' 
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Table  2.  Continued 


2«2X«*0EL  THTA**2Xf •PHI/FIT»i3X.»PHI/EXP*i3Xf 'DEL  PHI**/) 

100  F0PMAT(SXfF10.6f9F10.3) 

109  FORMAT  <lXt***»*  CONVERGENCE  CRITERIA  NOT  SATISFIED  AFTER'tlA* 
llXtMTEHATIONS  •••♦•) 

110  FORMAT<1H1*1x«*SUM  of  RESI0UALS*tE13.6*5Xf 'NEW  PROBABLE  ERROR*** 
1E13*6*5X**0L0  PROBABLE  ERROR  ■**E13*6) 

111  FORMAT(AFIO.A) 

112  FORMAT<1X*70(***) */*lX**««**««<»»  ERROR  *•*•*•*•**/* IX* 

l«NO  OF  CONSTANTS  ALLOWED  TO  VARY  EXCEEDS  TOTAL  NO  OF  DATA  POINTS 
2USEO**/«1X*70(***)) 

200  FOMMAT(1X*30(»»*).*  AERODYNAMIC  COEFFICIENT  EXTRACTION  BY  CHAPMAN 
IKIRK  method  *«30(***)«/  •1X**NSWTCH  s**I5**<  «0*  NORMAL  DATA  REDD 
2CTI0N*«/*17X**A1(  3D0F  DATA  GENERATED*  NO  ITERATION**/* 

317X*t«2«  3DOF  DATA  GENERATED  AND  STORED  ON  7S  PACK) • */» IX* * lUPOAT 
4***I5**(*0*  INITIAL  GUESSES  NOT  UPDATED  FOR  FOLLOWING  SHOT**/* 
SlTXftNE  0*  INITIAL  GUESSES  UPDATED)') 

302  format (IX** IaNGLE  »**I5«*(>0*  ALPA  AND  BETA  IN  TERMS  OF  VELOCITY  R 
1AT10S***/*17X**NE  0*  ALPA  AND  BETA  IN  TERMS  OF  ANGLE  OF  ATTACK  AND 

3 SIDESLIF|**/*1X**INLC0N  «**I5**(a  0*  INITIAL  CONDITIONS  EVALUATED 

4 FROM  OATA**/«17X**NE  0*INPUT  VALUES  USED)*) 

201  F0RMAT(  /*1X**PARAMETERS  ••••****»/«1X»*N0  of  points  (N)  ***15* 

I'TOTAL  no  variables  (NO  >**IS*5X**TOTAL  NO  ITERATIONS  (NITER)*** 
2I5*/«1X*'INTE6RATI0N  STEP  SIZE (DT) ** «EB.2*5X* *CONVERGENCE  BOUND  ON 
3 ITERATI'iN  (BOUND) *• *EB*2«//*1X* *BODY  ••••••*••*/« IX* 'AREA  *** 

4F10.6**FT**2**5X**OIA  * * tFl 0 .6* *FT * */* IX*  * IX  «* *F12.7* *SLUGS*FT**2 
5'*5X**IY  •*tFl2.7**SLUGS-Fr**2*«/*lX«*IZ  ■* *F12,7. *SLUSS-FT'**2' * 
65X«MZX  a'*F12*7**SLUQS-FT**2*) 

202  F0RMAT(/*1X**FREE  STREAM  *•***••***/* IX* *DENS1TY  (RO)  ■**E14.7* 
l'SLUbS/FT»*3**5X**U  ***F10.4* *FT/SEC* */  » iX* *COEFFICIENTS  ••*•••») 

203  F0RMATI1X*2(*AC0N(**I3**)  «* *E13.6*2X* * ( * *A8* * ) * *5X* *NC0N ( * * 13* 

1')  a**I2*5X)) 

301  format (IHI) 

303  format ( IX* *BEG1NNING  RECORD  NUMBER  «**1S) 

401  FORMAT(IOAS) 

402  FORMAT(1x*10A8) 

404  FORMAT(1X*40(***)**  INPUT  CARD  IMAGES  • *40( *•* I *//> 

406  FORMAT(lHi) 

407  FORMAT(lXt*NOTE  - **10A8) 

502  F0HMAT(1X*69(***) */*lX**'»*««  300F  GENERATED  DATA  STORED  BEGINNING 
20N  RECORD  NUMBER* * I5*2X* •••••**/* IX* 'UNDER  IGPNO  *i*I5*2X* 

3'  AND  IBlKNO  ***I5*/*1X**1I*1*T'*/*1X*'1I*2*PH1'*/*1X* 
4*II«3*THTA'*/*1X**1I*4*PSI**/*1X*69(***) ) 

WRITE(06*404) 

405  REaD(05*401<ENO*403)CARO 
WRITE(06*402)CARD 

GO  TO  405 

403  REWIND  5 
REAO(05*401)CARD 
WRITE(06*406) 

DEFINE  file  07(1140*4016«L* IREC) 

INITIALIZATION 

IPRlNTsO 

NPR1NT*0 
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Table  2.  Continued 


998  READ(0S«NAMtENDa999) 

AAC0N(6a>»IZX 

n*ntot 

WRlTE<06f301) 

WRITE  < 06*200) NSWTCHt lUPOAT 
WR1TE(06«302)  lANGLE* INLCON 
WRITE(06*A07)CARD 

WRITE t 06*201 )NTOT»NC «NITER*DTMIN>BOUND*AREAt 01 A* IX tlYtlZtlZX 

WR1TE(06*202)RO*U 

NC02>NC/2 

DO  204  I«1*NC02 

ITal^NCOZ 

AC0NI1T)«AAC0N(IT) 

ACONIDsAACONII) 

204  WRITE(06*203l I*ACON(l) *ALIST(I) *I*NCON(I) *IT*AC0N(1T>  * 
lALISTCIT) *1T*NC0N(IT) 

PI  " 3.i41592O«0 
A m PI/180. D«0 
1F(NSWTCH)8*8*7 

DATA  INPUT 

7 CONTINUE 
T<ll»0.D*0 
DO  10  laZfNTOT 
T<ll«T(l-l)*OTGEN 
PHlM(I)aU.D*0 
THTAM(I)*0.D«0 
10  PSlM(I)a0.D*0 
PHIM(l)a0.D«0 
THTAM(1>«0*0«0 
PSIH(l>a0.D«0 
GO  TO  9 

a continue 

IRECal 

500  REAO(07'lREC)16PT*lBtK«.II*TS 
IF(IGPNO.NE.IGPT)  GO  TO  500 
IFdBLK.NE.IBLKNO)  GO  TO  500 
IFdl.NE.l)  GO  TO  500 
REAO(07'lREC>lGPNOtIBLKNO«lI*PHIS 
HEAD (07 • IREC) IGPNOt IBLKNO* 1 I *THTAS 
REA0(07'1REC) IGPNO* IBLKNO* 11 *PSIS 

DATA  alteration 

DO  504  laltNTOT 
T(l)aTSdl 
PHIM(I)BPHlS(n 
THTAM(l)aTHTAS(l) 

504  P$IM<I)aPSIS(n 

lF(NPRINT.Ea. IGPNO. ANO.IPRINT.EO. IBLKNO)  60  TO  4 

WRITE<06*100)IGPNO«IBLKNO 

lRCTHPalREC-4 

WRITE (06*303)  IRCTMP 

LOOPaO 
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Table  2.  Continued 


NL00P«1 

DO  4 laltNTOT 

WRlTE(06tl02)T(I)*PSlH(I)  tTHTAMd)  tPHIMd) 

LOOP»LOOP*1 
KLO0PbNLU0P*50 
lF(L00P*Lr.KL00P)  60  TO  4 
NLOOPaNLUOP«l 

MRITE(06tlOO)  IGPNOtIBLKNO 
4 continue 

NPRlNTalGPNO 

IPRlNTalBLKNO 

call  part (SAMPLE tNDELTAfN«NTOT«NST ART) 

IF (NSESXH.CQ.l ) CALL  NOISE (SPHI «STHTA*SPSI (AHPHI «AMTHTAtAMPSI » 
lN*sAMPLEtNDELTA.NTOT) 

IFdNLCON.CQ.O)  CALL  ICS(ACON) 

00  46  I>ltN 
PSIMd)  a PSIMd)*A 
THTAMd)  a THTAMd)*A 
46  PHIMd)  a PHIHCnaA 
9 CONTINUE 
NSUMbO 

DO  11  I-ltNC 

11  IF(NCONd)  .EQ.l)  NSUMaNSUM«l 
SUM  a NSUH 

DEFINE  PARAMETERS  AND  TRANSFORM  COEFFICIENTS  IN  SUB  REG 

DO  12  Ial,NC 

12  CCONd)aNCONd) 

QINF  a RO*(U»»2)/2.0»0 

JXalZX/lX 

JYalZX/lY 

JZalZX/U 

lXPadZ-lY)/lX 

IYpadX-lZ)/lY 

IZP«dY-IX)/lZ 

SQTTaOIA/(2.D*0aU) 

SQT>QINFaAREA*DIA 

MXiaSOT/IX 

MX2«MXl*SQTT 

MYiaSOT/lY 

MY2aMYl*SQTT 

MZl-SOT/IZ 

MZ2»MZl*5aTT 

call  RE6(AC0N) 

XNaN 

CHECKa3.U»0«XN-SUM 

IF (CHECK. LE. 0.0*0)  60  TO  997 

NTaNC*! 

EaO.O*0 

NRaO 

20  continue 

JXaCUN(6U)/lX 

JYaC0N(6U)/IY 

JZaC0N(6U)/IZ 
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integration 

call  INTQRL(FCT*ORKGStNtERROR»AUXfY*DERY»DTMIN) 
IF(NSWTCM,6T.0>  60  TO  *7 
IF<NSUH*EQtOI  GO  TO  47 
NR«NR*1 

CHECK  FOR  CONVERGENCE 

EUOSQRT(RESSUM/(3.0«0*XN-SUMt ) 

RESSUHaOSQRT (RESSUM) 

MRlTE(06«110)RESSUMtEl*E 
IF (DABS (El-E> -BOUND) 36t36t 21 

21  IF(NR-NlTER)22«22t36 

22  NTRiGaO 
38  continue 

IFlNSUHtEO.l)  GO  TO  33 
E>E1 

contract  C matrix  • REMOVE  ZERO  ROWS  AND  COLUMNS 

DO  23  I«1»NC 
DO  23  jaltNT 

23  CTdf J)«C(ltJ) 

LL*0 

DO  90  lalvNC 

92  IF!  <I«LL) •GT.NC)  60  TO  90 
IF|NC0N(l*LL)>91t91«93 

91  LLbLL«1 
60  TO  92 

93  L«LL*I 
DO  94  K«ltNT 

94  C(I*K)aCT(LtK) 

90  continue 

LL«0 
LMaNSUM 
DO  80  I*ltNC 

82  1F( (I*LL) •GT.NC)  GO  TO  80 
IF(NC0N(l*LL>>61*81tB3 
81  LL«LL*1 
GO  TO  82 
B3  LaLL^I 

DO  84  K«1»LM 

84  CT(K«I>«C(KfL) 

80  CONTINUE 

DO  85  I«ltLM 

85  CT(l*LM*'l)«C<ltNT) 

L«NSUM 
LLbNSUM*! 

DO  71  I«ltL  ^ 

DO  7l  jaltL 
71  C(lf J)aCT(ltJ) 

C 

c matrix  inversion  of  c 
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IF(NTRIQ)40t27»40 

27  call  INV(CT«LtLLtWX*CC> 

GO  TO  A1 

AO  CALL  INV(CT«L*L*MX«CC) 

A1  DO  32  l^ltNC 
C(1»LLI«0.0«0 
WX(l)«CT(Iil) 

32  MXX<I)«0tO«0 
LHaO 

EXPAND  C INVERSE  TO  INCLUDE  ZERO  ROWS  AND  COLUMNS 

DO  31  I«1*NC 
lF(NCONU))28t2Sf30 

28  LMaLH-1 
GO  TO  31 

30  ITbLH«I 
CntLL)«CT(lTtLL> 

WXX(I)aWX(IT) 

31  continue 
GO  TO  34 

SOLVE  FOR  CORRECTION  TO  COEFFICIENTS 

33  00  35  I-ltNC 
IF|NC0N(I)-1I39«49*49 

49  C<I»NT)«C(UNT)/C(ItII 
NXX(I)al.0«0/C(ItI) 

60  TO  35 

39  C(I*NT>«0.0«0 
NXX(l>a0*D«0 

35  continue 
LL*NT 

34  1F(NTRIG.EQ.1I  GO  TO  42 
WRITC(06tl04)NR 

DO  60  I*ltNC 
' 60  WXX<1I«C(I*LL) 

TRANSFORM  COEFFICIENTS  AND  CORRECTIONS  IN  SUB  REGOUT 

CALL  REGOUT (ACONtWXX) 

DO  50  I«1»NC02 
ITal«NC02 

50  WRlTC(06tl03)ALIST(I> tACON(I) tWXX(I) fALISTlIT) f ACON ( IT) t NXX ( IT) 
00  43  laltNC 

43  CON(I)aCON(I)«C<l»LL) 

60  TO  20 

36  CONTINUE 
NTRlGal 

37  60  TO  38 
42  continue 

cumpute  probable  error  in  coefficients 
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DO  44  l«ltNC 

44  tfXX(n«El*OSGRT(OABS<WXX(I) ) ) 
call  RC60UT<ACONtWXXI 
WRlTE(06tl05) 

1F(NR.GT  .NlTER>WRITE(06tl09)NR 

DO  45  I«lfNC02 

lTsl«NC02 

45  MRITE(06tl06)  1« ACON ( I ) tWXX ( 1 ) t ALIST (I > t ITt AC0N( IT) «WXX < IT) 
IfALlSTdT) 

47  WRlTE(06fl07) 

LOOPsO 

NLOOPal 

IRECalRCGEN 

DTTEMP«T(2)-T(1) 

DO  46  I>1»N 
TT»T(1) 

21  " PSKD/A 

22  * PSIM(I)/A 

23  > THTA(I)/A 
Z4  ■ THTAMiD/A 

25  > PHI  (D/A 

26  « PH1M(1)/A 
2lO«Zl-22 
ZD-Z3-Z4 
Z12-Z5-Z6 

WRlTE(06«l08)TTtZl«Z2tZ10«Z3iZ4fZll«Z5tZ6tZ12 

THTAS(I)«Z3 

PS1S(D>Z1 

TS(I)aTT 

PHIS(I)«2S 

IF(NSWTCH,NE.O)  GO  TO  597 

TS(1)>Z10 

PHIS(I)>211 

597  continue 
LOOP«LOOP*l 
KLOOPsNLOOPASO 
IF(LOOP.LT.KLOOP)  GO  TO  46 
NLOOPaNLOOP*! 

WRITE(06tl07) 

48  CONTINUE 

IFiNSWTCH.NE.O)  GO  TO  599 

COMPARISON  OF  FIT  DATA  AND  EXPERIMENT  STORED 
IRECalRCSTR 

lF(lRCSTR.6EtIRC3EN-3.AND.IRCSTR.LE.IRCGEN»3)  60  TO  599 
WRITE (06t598)NSMTCHtlRCSTR«I6PNO«IBLKNOtOTTEMP 

598  FORMATdHltlXt'NSWTCH  * I5i/t 1X« 'SOOF  DATA  STORED  BEGINNING  AT  RE 
ICORD  t,I:},2Xt*UNDER  IGPNO  >• tI5t2X« 'AND  IBLKNO  • I5»/* IXt • I I»5tTH 
2TA»i/tlX»»II«6tTHTA-THTAM»t/tlXt*lI*7*PSI'i/»lX»Ml*8tPSI-PSIMN 
3/tIX**DT  ■•*F7*4i2X*'SEC>i1H1) 

llaS 

WRITE(07'1REC) IGPNO* IBLKNOt I I «THTAS 
lla6 

WRITE (07<IREC) IGPNO* IBLKNO* I I tPHlS 
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n«7 

wRirE(07*IREC)16PNOtlBLKNOtIl,PSIS 

n>8 

WRITE ( 07 'IREO IGPNOtIBLKNOi II*TS 
599  continue 

IF(NSNTCn.NE.2)  GO  TO  501 

3DOF  GENERATED  SOLUTION  STORED 

MRl TE ( 06  tS02) IRCGENt IGPNOt IBLKNO 
1I«1 

WRITE ( 07 'IREO IGPNOtIBLKNOtlItTS 
ll»Z 

WRITE <07 'IREO IGPNOtIBLKNOtlltPHIS 
II>3 

WRITE  1 07  < IREC) IGPNO* IBLKNO* 1 1 tTHTAS 
II-A 

WRITE <07 'IREC) IGPNO* IBLKNO* I I tPSlS 
501  continue 
GO  TO  996 
997  WRITE<06*112) 

996  continue 

IF<IUPOAT.EO.O>  GO  TO  996 
DO  995  I>1*NC 
995  AACON<I>aACON<I) 

IZX«AAC0N<60) 

60  TO  998 
999  CONTINUE 
STOP 
ENO 


SUBROUTINE  NOISE <SPH1*STHTA*SPSI*AMPHI*AMTHTA*AMPS1*N* 
1SAMPLE*NDELTA*NT0T) 

THIS  SUBROUTINE  ADOS  GAUSSIAN  NOISE  TO  DATA 
IMPLICIT  REAL*8<A-H*0*Z> 

COMMON/MEAS/PHIM<1000) *THTAM<1000) «PSIH<1000I 
COMMON/S0L/T<1000) *PHI <1000) .THTA < 1 000 ) (PSI < 1 000) 
WRITE<06*100>SPS1*STHTA*SPHI«AMPSI*AMTHTA*AMPHI*N 
call  TIME<IYEAR*ISEC) 

7 1F<1SEC-1000)  5*5*6 
6 1SEC>ISEC-1000 
GO  TO  7 

5 lXa<lSEC«lYEAR)*2*l 
lYalx*2 
IZ>lX*4 
S1G1PH«0*0«0 
SlG2PHaO.O«0 
SI63PHa0.O«O 
SlGlTHaO.O*0 
Sl62THaO.O«0 
SlG3THaO.O«0 
S16lPSaO.O*0 
SlGZPSsO.O^O 
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S1G3PS«0«0«0 

TOTPHl-O.O^O 

tottha«o.o*o 

TOTPSI»0.O«0 

LOOPal 

NLOOP«l 

00  10  I«1*N 

CALL  OAUSSdXtSPHI.AMPHl.VPHI) 
call  gauss ( lYiSTHTAt AMTHTA. VTHTA) 
call  GAUSS(lZ*SPSl*AMPSltVPSl) 

PHIM(l)aPHIM(l)*VPHl 
THTAM ( 1 ) »THTAM 1 1 > ♦VTHTA 
PS1H<I)«PSIM<1)4VPS1 
TOTPHI«TOTPHI^VPHI 

tottha-tottha*vthta 

TOTPSI»TOTPSI*VPSI 

P«DABS(VPHn 

0>DABS(VTHTA) 

RsDABS(VPSI> 

lF(P.LTiSPHl)SIGlPH>SIGlPH^UD«0 

IF(P.6E.SPH1.ANO«P.LT.2«0*0«SPH1>SI62PHbSIG2PH+1.0^0 

1F(P.GE.2.0*0*SPHI,AND.P.LE,3.040*SPHHSIG3PH«SIG3PH^I,0^0 

IF(Q.LT.STHTAISIGITH»S1GITH^I,040 

IF{0.6E»STHTA.AND.Q.LT*2.0^0*STHTA)SI62TH»SIG2TH+I.D^O 

lF<C».6E.2,0^0»STHTA.ANO.a.LE.3.D*0«STHTA)SI03TH-SlG3TH4l.O^O 


lF(R.LTtSPSnSlGlPS«SIGlPS4l.D40 

lF{R.6EtSPSI.AN0.R.LT.2.O40*SPSI>S162PS»SIG2PS*l.D40 
IF(R.GE.2.040»SPSI.ANO.R.LE.3.040»SPSI)SI63PS»SIG3PS4l.D40 
10  CONTINUE 

WRITE (06tl03>TOTPSItTOTTHA»TOTPHItSlGlPS»SI61TH»SIOlPHf 
lSIG2PS«SlG2THfSlG2PH*SlG3PStSlG3THtSIG3PH 
100  FORMATdHlflX*'***  GAUSSIAN  NOISE  ADDED  ••••*/»lXt 

PONE  STANDARD  DEVIATION  - PSI  «« tFlO.GtSXi 'THETA  ■•fFlOtStSXi 
2'PHI  ■•«F10.6»/flXt'CONSTANT  BIAS  - PSI  •NFlOtOfSXi 
3'THETA  ■•*Fl0.6f5X»»PHI  ■••Fl0.6*/flXf 'TOTAL  NO  OF  TIME  POINTS  USE 
AD  IN  FIT  (N)«'tI5) 

103  format IlHltlXt 'SUM  OF  NOISE  ADDED* t/fbXf *PSI ' »6X» 'THTA* iTXt 'PHI ' t 
l.tlXt3F10,3*//flXt'NOISE  CONTRIBUTION  WITHIN  ONE  SiGMA'#/* 
26X»'PSI' t6Xi •THTA'»7X»*PHI'f/»lXt3F10.3i//f 
31X« 'NOISE  CONTRIBUTION  BETWEEN  ONE  SIGMA  AND  TWO  SIGMA'*/* 

♦6X* 'PSI **6X.*THTA**7X. 'PHI *,/«lX*3F10.3*//»lX* 'NOISE  CONTRIBUTION 
5BETWEEN  TWO  SIGMA  AND  THREE  SIGMA' */f6X* *PSI • t6X* *THTA* tTX* *PHI * * 
6/tlX*3F10,3) 

RETURN 

END 


SUBROUTINE  RANDUlIXtlYtTPL) 

IMPLICIT  REAL*8(A-H*0-Z) 

IY»IX*65539 

IF(lY>5i6*6 

5 IYbIY*21A7A836A7*1 

6 YPL«IY 

YPL*YPL«.A6566130-9 
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RETURN 

END 


SUBROUTINE  GAUSS! 1X«S«AM*V) 

IMPLICIT  REAL*e(A-HtO-Z) 

Aa0«0*0 

00  SO  I>I*A8 

CALL  RANOUdXvIYtV) 

IX»IV 
50  A»A*Y 

V«<A-2A,U*0)/2.0^0 

V»V*S*AM 

RETURN 

END 


SUBROUTINE  ICSIACON) 

IMPLICIT  REAL*8  (A-H*0-Z) 

DIMENSION  AC0NI72) 

C0MM0N/MEAS/PHIM(1000)  iTHTAMdOOO)  iPSIMdODO 
COMMON/SOL/T (1000) tPHI (1000) tTHTA d 000) «PSI dOOO) 
UT«T(2)-Td) 

PHIO-PHlMd) 

THTAO>THTAMd) 

PSIOsPSlMd) 

PHIP  ■ (PHIM(2)-PHI0)/DT 
THTAP  ■ (THTAM(2)-THTA0)/0T 
PSIP»(PS1M(2)-PSI0)/DT 
ACON(70)>PH1P«PS1P*DS1N(THTAO) 

ACON(71)sTHTAP*OCOS(PHlO)«PSIP*DCOS(THTAO»*OSIN(PH10) 

ACON (72) "PSIP^OCOS (THTAO) *0005 (PHIO) -THTAP*0SIN (PHIO) 

ACON(67)aPHIO 

AC0N(68>>THTA0 

AC0N(69)aPSI0 

RETURN 

END 


SUBROUTINE  PART (SAMPLEtNOELTA*NtNTOT«NSTART) 

C 

IMPLICIT  REAL*8  (A-H*0-Z) 

COMMUN/MEAS/PHIMCIOOO)  tTHTAMdOOOl  (PSIMCIOOO) 
COMMON/SOL/T (1000) tPHl ( 1 000 ) t THTA ( 1 000 ) *PSI d 000 ) 
DO  13  l>2tNTOT 
1F(T(1) .NE«0.0*0)  NTEMPal 
13  continue 

TMAXbT(NTEMP) 

DO  11  I>NSTART»NTOT 
J»d-1)«NDELTA*1 
IF(J.GT.NTEMP)  GO  TO  11 
RATIO* ( T ( J » -T (NSTART) ) /TMAX 
IF (Ratio. GT. SAMPLE)  60  TO  11 
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N«I 

T(1>bTIJ) 

PHlM(I)aPHlM(J) 

thtam<I)»thtamij> 

PSIM(I)aPSlH(J) 

11  continue 

NbN^I-NSTART 


DO  12  1>1«N 
J»I*NSTAHT-l 
T(I)aT(J> 
PMIM(n»PHIM<J» 
THTAM(I>«THTAM(J» 
12  PSIM<I)aPSIM(J) 


15 

lA 

10 


N1>N«1 

IF(NTOTtLT.Nl>  GO  TO  lA 
DO  15  laNltNTOT 
TCI)«0.0*0 

continue 

WRlTE(06tl0)SAMPLEtNDELTAtNSTARTtNfT(l>  tT (N) 
FORMAT(lHl»lXt*«*»*  DATA  BLOCK  ALTERED  •**•» t/f 1X» 'SAMPLE  ■' 
lF7,A,/f IXi 'NDELTA  »itl5*/»lX*'NSTART  ■*tI5i/tlXt'N 
21Xi*T(ll  ■ •*F7*A»SXt*T(N>  ■ *tF7.A*/» 

3lXf28(»»M) 

RETURN 

END 


SUBROUTINE  ORKGS<XI*VfDERYiNOIMtFCT«AUXtOT*NOT) 

IMPLICIT  REAL»8U-H»0-Z) 

DIMENSION  Y(l) iDERYCl) lAUXt  5*1) 

XaXl 

DO  5 JalvNOT 
DO  1 laIvNOlM 

1 AUX(5*I)«Y(1) 
call  FCT(X*Y.0ERY) 

DO  2 lalfNOlM 
AUX(1*I>«DERYII) 

2 Yd)  >AUX  <5*  1 ) « •SD4'0*DT*AUX  < I « I ) 

XTEMPbX*.5O*0*OT 

call  FCT(XTEMP*Y.0ERY) 

DO  3 Ial*NOIM 
AUX(2tI)»0ERY(l) 

3 Y (I ) >AUX (5* 1 ) « .SO*0*DT*AUX (2* I ) 

CALL  FCT(XTCHP*Y*DERY) 

DO  A lalvNOlM 
AUX(3«I)aOERYd) 

A Yd)«AUX(5*I)'>DT*AUX(3*l) 

XaX*DT 

call  FCT(X»Y.0ERY) 

DO  5 laltNDIH 

AUX(A*I)aDERY(I)  , ^ , 

5 Yd)«AOXJ5il)*DT»«AUX(l*l)*(AUX<2»l)*AUX(3*I>>«2.D-»0 
l«AUX<Ail) )/6.0«0 
RETURN 
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END 

c 

C 

FUNCTION  OUHIltJ) 

C 

IMPLICIT  REAL*6(A-HtO-Z> 
DUM*0.D*O 

IFd.EQtU)  DUMsl«D»0 

RETURN 

END 


SUBROUTINE  INV (Cf NCtNCSl tMXXtCC) 

IMPLICIT  REAL*e(A>HtO-Z) 

DIMENSION  C(72f 731 tWXX (72) «PIVOT (2) «CC (72tl44) 
3000  FORMAT(/10X«*DET  IS  EQUAL  TO  ZERO*) 

NCT*NC*2 
NCP1>NC*1 
DO  10  I«1*NC 
DO  lO  JaltNC 
10  CCdf J)>C(I*J) 

DO  20  I>itNC 
DO  20  jaNCPltNCT 
20  CCdtJ)«0,0*0 
DO  30  laltNC 
30  CCd«NC*l)«l.D^O 
DO  205  laltNC 
PlVOTID^CCdtl) 

DO  200  K-ltNC 
PlV0T(2)aCC(Kf 1) 

126  IF(K.I)  13Stl30tl40 
130  DO  150  JaltNCT 

IFCPIVOTCD)  134*210tl34 

134  CC(KtJ)>CC(IiJ>/PlVOTdl 
150  continue 

60  TO  200 

135  DO  160  JaltNCT 
IFCPlVOTd))  136tl60«136 

136  CC(KtJ)>CC(K«J>*CC(If J)*P1V0T(2)/PIV0T(1) 

160  CONTINUE 

GO  TO  200 
140  DO  170  JaltNCT 

IF(PIV0T<Z))  145fl70d4S 
145  CC(K,j)aCC(Kt J)/PIVOTI2)-CC(I*J) 

170  CONTINUE 
200  CONTINUE 
205  CONTINUE 
GO  TO  250 
210  WRITE  (6i3000> 

250  DO  300  lalfNC 
DO  300  Jsl*NC 
300  C(ItJ)sCCdtJ«NC) 

NCSaNCSl 

350  IF(NCS-NC)500«500t400 
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400  DO  420  1>1*NC 
MXX(l)>C(ItNCS) 

420  C(ItNCS)s0.D*0 
00  450  I>1*NC 
DO  450  JaltNC 

4S0  C( 1 *NCS)»C I ItNCSl ♦€ (It J)*WXX IJ) 
NCS«NCS-1 
GO  TO  350 
500  continue 
RETURN 
END 


SUBROUTINE  1NT6RL (FCTi0RKGS«N«ERR0R*AUX*Y»0ERYi0TMIN> 

IMPLICIT  REAL*e(A-HiO-Z) 

DIMENSION  AUX(  5 t 1 > t Y ( 1 ) tDERY ( 1 > 

COMMON/MEAS/PHIM(1000)  tTHTAMIlOOOl  tPSIMdOOO) 

COMMON/SOL/T(lOOO»  tPHl  (1000)  tTHTAdOOO)  «PSI  dOOO) 

C0MM0N/C0EF/C0N(72) fCC0N(72) tNCiNC0N(72l 

C0MM0N/MATRIX/C(72*73>  tRESSUMiSPHI ♦STHAtSPSl 

Aa6*283iaS310«0 

NLsNC-5 

00  10  Ialt6 

10  Y(I>aC0N(NL«I-l) 

DO  11  laltNC 
Y(6*I)aOUM(NL«l) 

DO  11  JaliS 

11  Y(6*J*NC«l>aOUM(NL«J«I) 

NEQa(NC*l)*6 

ONEOaNEQ 

NLaNC«l 

Kal 

PHl(K>aYd) 

THTA(K)aY(2) 

PSI<K)sV(3) 

THTA(K)«THTA(K)«C0N(6SI 

PSl(K)aPSI<K)«C0N(66) 

OPHlaPHlM(K)-PHl(K) 

lF(OPHI.LTt(-A/2.D*0))  PHI (K) aPMl (K) -A 
IF(OPHI.GT.(A/2.D*0)>  PHI (K» aPHI (K) ♦A 
DO  19  laliNC 

Y(6«NC»I)aY(6«NC*ll«0UH(65iI> 

Y(6*2*NC«I)aY(6«2»NC«I)«0UM(66fI) 

19  CONTINUE 
DO  12  lalfNC 
DO  12  jaltNC 

C(ltJ)aY(6*I)*Y(6*J)/SPHI*a2  ♦ Y(6*NC*I) •Y(6*NC*J>/STHA**2 
!♦  Y (6*2aNC*I I ^Y (6*2«NC*J) /SPS1»»Z 

12  Cd*J)aCd»J)aCCONd>*CCON<J) 

00  17  jaltNC 

C(J»NL)a(PHIM(K)-PHI (K) ) ^Y (6* J) /SPHl**2  ♦ (THTAM (K) -THTA (K) ) • 
1Y(6*NC*U)/STHA**2  ♦ (PSIM (K) -PSI (K) ) *Y (6*2aNC« J) /SPS1**2 
17  C(J*NL)aC(JtNL)*CCON(J) 

RESSUM  a ((PHlM(K)-PHI(K))/SPHl)*a2  ♦ ( (THTAM(K)-THTA(K) )/STHA)**2 
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1*  ((PSIM(K)-PSl(K))/SPSn**2 
00  13  KBZfN 
XlaT(K<-l) 

DTEBT(K>-T(K-n 

OTaOTE 

NOT-1 

lF(UTHIN.GE.DTei  GO  TO  18 

N0T»1)TE/0TMIN 

OTaNOT 

OTaOTE/OT 

18  continue 

call  ORKGS(XItYtOERY«NEQ«FCTfAUXiOTtNOT) 

IF<YI1) .LE.(«A> )Y(l)aY(l)«A 
IF<Y(1) .GE.AIY(l>aY(l)oA 
PHI  (K)aY(n 
THTA(K)ar(2) 

PSI (K)aY(3) 

THTA(K)aTHrA(K) »C0N<65) 

PSI (K)aPSI (K)«C0N(66I 
DPHlapHIM(K)<-PHI  (K) 

IFtOPHI.LT. (-A/2tO*0> ) PHI <K|aPHI (K)-A 
1F(0PH1.6T. (A/2. 0*0) ) PHI (K)aPHI (K) *A 
00  20  laltNC 

Y (6*NC*I ) aY (6*NC*I ) ♦OUM (65* I ) 

Y (6*2BNC* I ) ay (6*2«NC* I > *0UH (66. I > 

20  CONTINUE 
DO  15  lal.NC 
00  15  jal.NC 

C(I.J)aT(6*l)aY(6*J)/SPHl**2  • Y (6*NC*I > ^Y (6*NC*J) /STHAaa2 
1*  Y(6*2*NC*1)»Y(6*2*NC*J)/SPS1**2  * C(I.J> 

15  Cdf J)aC(ItJ)*CCON(I)*CCON(J> 

DO  16  jaliNC 

C(J.NL)a(PHlM(K)-PHl(K) )*Y(6*J)/SPHl*a2  ♦ (THTAM (K ) -THTA (K) ) • 
1Y(6*NC*J)/STHA**2  ♦ (PSlH(K ) -PSI (K) ) *Y (6*2»NC* J> /SPSI*»2  ♦ C(JtNL) 

16  C(J*NL)«C(JtNL)*CCON(J) 

RESSUM  a ( (PHIH(K)-PHI (K) )/SPHI)**2  * ( (THTAM(K>-THTA(K) )/STHA) *#2 
1*  ((PS1M(K)-PSI(K)>/SPS1)**2  * RESSUH 
13  continue 
RETURN 
END 
C 
C 

SUBROUTINE  MBARA (HI *M2.H3iPSI .THTA.PHI tP.Q.R.Pl tP2.P3t 
1 Q1.Q2.03.J) 

C 

IMPLICIT  REALMS (A-H.O-Z) 

KEAL«8  MAGfMl«M2.M3«MYPX.MYPZtMYTP*MZPXtMZPYiMZTPt 
1 MYT.MZT.MX.MY.HYZ.MZfMZYtMXP.MYP.MZP.HYHZ.HZHY 
C0MH0N/TC0EF/A.B.MA6*ALPA»BTA«DLTA2.RHAt0MAtRHBtDMB«RHABtDMAB 
C0MM0N/C0EF/C0N(72) fCC0N(72) «NCtNC0N(72) 

C0MM0N/C0NST2/RX1 .RX2.RY 1 tRY2«RZl .RZ2. AA 

TEMPYlaCON(17>*CON(18)*ALPA*CON(19)*RMA*CON(20>*ALPA**2*CON(2l)*AL 
1PA*RMA*C0N(22)*RMB*C0N(23)*BTA**2*C0N(2A)*ALPA*RMB*C0N(25)*RMAB 
TEMPYZaOUH ( 17 . J>  *DUM ( 18. J> •ALPA*CON (18) •A*OUM ( 19. J) •RMA*C0N ( 1 9) aOM 
lA  «OUM(20.J)*ALPA**2*2.D*0aCON(20)aALPA*A*OUM(21.J)*ALPA*RMA 
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♦CON  121 ) • < ALPA»DMA+RMA»A) ♦DUN  <22f J>  *RMB^CON (22) *DMB 
♦OUH (23 t J) ♦BTA*»2^2tD^O*CON ( 23) *BTA»B^0UM (24 tj) •ALPA*RMB 
♦CUN(24)* (ALPA»OMB^RMB»A)  ♦ DUM(25« J) »RMAB^C0N(25} •DMAS 
«DUM(ltJ)^DUM(2*J)«ALPA^C0N(2)*A  ♦ DUMO* J)»RMA^CON(3)«OMA 
♦DUM (4t J) *ALPA««2^2*0^0*CON (4) •ALPA*A^DUM (5t J) •ALPA*RMA 
♦CON (S) * ( ALPA*OHA^RMA*A) ♦DUM (6* J) •AUPA**3^3.0^0*CON(6) *4 
•AUPA*«2  ♦DUM ( 7 * J) •RMA*ALPA*«2^C0N (7 ) • (0HA»ALPA**2 
♦2«0^0*A*ALPA*RMA)«OUM(6«J)*RMB  ♦ CON (S) mONB^OUM (9* J) • 
BTA**2^  2.D^0«CON(9)*BTA*B^DUN(10»J)*RMB*BTA**2 
♦C0N(10)*(OMB*BTA*»2  ♦ 2.D«0*BTA*RHB*B) ♦OUM ( 1 1 • J) •ALPA*RMB 
M2«TEMPY^  CON(ll)*(ALPA*OMB^A*RMB)  ♦OUM ( 12t J) •RMAB^CON ( 12) «OHAB 

1 ♦DUM (13*0) •RMB*ALPA*»2^C0N ( 13) * (DMB»ALPA*»2^2.D^0»ALPA*A*RMB) 

2 ♦DUH ( 14  * J) »RHAB»ALPA^CON ( 14 ) » ( ALPA*OHAB^RHAB*A) ♦DUM ( IS* J) *ALPA 

3 •BTA**2^C0N ( 15) • (2t0^0»ALPA*BTA*B^A*BTA**2) ♦DUM( I6t J) •RMA»BTA** 


42  ♦CON(ie)»(DMA*BTA»<»2^2.D^O*B7A*RMA*B)^  Q2*TeMPYl^Q*TEMPY2 

5 ♦O3*C0N(26)*eTA^R  •DUM(26fJI •BTA^R*C0N(26)*B  ♦Q1«C0N (27) *674 

6 *P«DUM(27«J> •BTA^P«C0N(27)*B 

TEMPZ 1 ■ CON (36) ♦CON ( 37 ) •ALPA^CON (38) •RMA^CON (39) *ALPA**2^C0N (40 ) * 
1RM8  ♦CON(41)*8TA**2^CON(42)*ALPA*RMA^CON(43)*ALPA«RMB*CON(44)»RMA8 
TEMP22«0UM(36*J)«DUM(37tJ)*ALPA^C0N(37)*A^0UM(38tJ>»RHA^C0N(3a)»0M 
lA  ♦OUH(39«J)*ALPA«*2*2*0^0*CON(39)*ALPA*A  ♦ OUM(40t J)*RMB 

2 ♦C0N(40)»0MB^DUM(4l»J)«BTA*»2  ♦ 2.D^O*CON (41 ) *BTA*B 

3 ♦DUM(42»J)*ALPA*RMA^C0N(42)*(ALPA«DMA^RMA*A)^DUH(43tJ) 

4 •ALPA*RHB^C0N (43) • ( ALPA*0HB^RMB*A) ♦DUH (44 1 J) •RMAB 

5 ♦C0N(44)*0MAB 

TEMPZbDUM (28* J) ♦DUM (29* J) •BTA^CON ( 29) ttB^DUM (30* J) *ALPA*BTA 

1 ♦CON(30)*(ALPA«B^BTA*A)^OUM(31*J)*RMA*BTA^CON(31)*(BrA*DMA 

2 ♦RMA*B) ♦OUH(32*J)*BTA**3^3.D^O«CON(32)«B*BTA**2 

3 ♦0UM(33*J)«BTA*ALPA**2  ♦ 2.D^O*CON(33)*ALPA«BTA*A^CON(33)*  B 

4 •ALPA^^Z  ♦ DUM(34*J)*BTA*RMB*C0N(34)*(BTA*0MB^RM8*B) 

5 *OUM(35iJ)*BTA*RMAB^CON(35)*(8TA«OMAB^RHAB*B) 

M3n  TEMPZ^  a3*TEMPZl ♦R*TEMPZ2^Q«BTA*DUM(45* J) ♦CON (45) • (02»BTA 
1 ♦Q*BI  ♦ DUM(46*J)«P*AUPA^C0N(46)»(Q1*ALPA^P*A) 

MIb  DUM(47*J)^  0UM(48*J>»ALPA^C0N(48)*A  ♦DUM (49*J) •BTA^CON (49) «B 

1 ♦DUM (50* J) •ALPA*BTA^C0N (50) • ( A»6TA^ALPA*B) ♦DUM (51 *J) •RMA*BTA 

2 ♦C0N(S1>*(BTA«0HA^RMA«B)  ♦0UH(S2(J)*BTA*RMB^C0N(52) 

3 •(BTA«0MB^RMB«B)  ♦ DUM(53*J)*P^C0N(53)*Q1 
THTAPbQ«DCOS (PHI ) -R«DS1N (PHI ) 

PSIP>  IR*OCOS (PHI ) ♦Q«0SIN (PHI) ) /DCOS (THTA) 

THTAPAaQ2*DCOS (PHI ) -Q*P1*DSIN(PHI ) 

1 >Q3*DSIN(PHI)-R*P1«0C0S(PHI) 

PSIPA»PSIP*DTAN (ThTA) •P2* (Q3*DC0S (PHI ) -R*P1«0SIN (PHI ) ♦ 

1 02*0SIN (PHI )^Q*P1«DC0S (PHI) ) /DCOS (THTA) 

MYTPaOUM (63* J) •THTA^CON (63) «P2 

MYTP»MYTP^0UM(62*JI»THTAP^C0N(62)»THTAPA 

MZTPbOUM (64* J) *PSI^CON (64) *P3 

MZTP«MZTP^DUH(61*J)*PSIP^C0N(61)*PSIPA 

MYT*C0N(63)«THTA 

MYT"MYT*C0N(62)«THTAP 

MZT"C0N(64)«PSI 

MZT>MZT^CON (61 ) •PSIP 

MYPA«MYTP*RX1/RYI 

MYPZ»MYTP«RZ1/RYI 

MZPXaM2TP*RXl/RZl 

MZPYbMZTP»RY1/RZI 
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MX«MYT*0S1N (PSD *0005 (THTA)*RX1/RY1-MZT*0SIN<THTA>*RX1/RZ1 

MY«MYT»0C0S(PSD 

MYZ"MY*RZl/RYl 

MZ*MYT*DSIN (PSI ) •OSIN (THTA) •RZl/RYl ♦MZT»0C0S( THTA) 
MZY*MYT*OSlN(PSD*DSlN(THTA)«HZT*DCOS(THrA)*RYl/RZl 
MXP»MYPX*DSIN (PSD  *0C0S (THTA) ♦MYT* (P3«DC0S (PSI ) •OCOS (THTA) - 
1 P2*I}SIN(PSI)*0S1N(THTA) )*RX1/RY1-MZPX*0SIN(THTA) 

2.MZ  T*P2«UCOS ( THT A ) •RX 1/RZ 1 
MYP»MYTP«0C0S (PSI ) -MYT*P3*0SIN (PSD 
MZP>MYPZ*0SIN (PSI ) •OSIN ( THT A) « 

1 MYT^(P3*OCOS(PS1)^OS1N(THTA)*P2^0SIN(PSD^OCOS(THTA)  ) 

2 •HZI/RY1*MZTP^0C0S(THTA)-MZT^P2^0SIN(THTA) 

M YHZ>MYPZ^0C0S (PSI) -MYT*P3*0S 1 N ( PSI ) •RZ 1 /R Y I 
MZHY»MYTP^0SIN(PSD^0SIN(THTA)* 

1 MYT^(P3^0COS(PS1)^OSIN(THTA)*P2^0SIN(PS1)^OCOS(THTA) ) 

2 «M2PY^0C0S( THTA) -MZT^PZ^OSIN (THTA) •RYl/RZl 
MlaMDMXP 

H2»M2-P 1 *M Y«OS IN ( PHI ) YP^OCOS ( PH I ) 

1 •MZY*P1^DC0S (PHI ) «MZHY*0S1N (PHI ) 
M3aM3«HYZ^PM0C0S(PHl)-MYHZ^DSIN(PHI) 

1 -MZ*P1^0SIN(PH1)«MZP^OCOS(PH1) 

M1bM1«q1^(C0N(S4)*RMA*C0N(5S)^RMB«C0N(56)^AUPA^^2«C0N(57)«BTA^^2) 

1 (RMA^OU)(  (S4t J) ♦OMA*CON (54) *RMB^OUM (55« J) «OMB^CON (SS) 

2 «OUM (56 t J) •ALPA^^2*2.0*O^CON (56 ) •ALPACA 

3 «OUM(57tJ)«BTA«^2*2tO«O^CON(57)*BTA^B) 

RETURN 

END 

C 

C 

SUBROUTINE  MBAR (MXOI *MY01 tHZOl iPSI *THTA«PH1 tP»0«R) 

C 

IMPLICIT  REALMS (A-HfO«Z) 

REALMS  MXOI tMYOI tHZOl tMAOtMYT (MZT 

COMMUN/TCOEF/AtBtHA6tALPAtBTAtOLTA2iRMAtOMAtRHBfOMBtRMABtOMAB 
C0MM0N/C0EF/C0N(72) tCC0N(72) tNCtNC0N(72) 
C0HH0N/C0NST2/RXltRX2tRYItRY2tRZltRZ2tAA 
MXOlaCON (47) tCON (48) •ALPA«CON (49) •BTA 

♦CON (SO) •ALPA^BTA*CON (51 ) •RMA^BTA 
♦CON (52) •BTA^RMB^CON (53) •P 

♦P^ (CON (54) aRMA^CON (55) •RMB^CON (56) •ALPA*^2^C0N (57) •BTA»*2) 
TEMPYaCON ( 17) ♦ CON ( 18) •ALPA^CON ( 19) •RMA^CON (20) *ALPA^^2^C0N (21 ) • 
ALPA*RMA^C0N (22) •RMB^CON (23) •BTA«*2^C0N (24) aALPA^RMB^CON (25) a 

rhab 

MYOl  a CON(l)aCON(2)aALPAaCON(3>aRMAaCON(4)aAl.PAaa2«cON(5)aALPA 
aHMAaC0N(6) aALPA^^3aC0N (7) aRMAaALPAaaZ^CON (8) aRMBaCON (9) 
aeTAaa2*cON(10)aRMB^eTAaa2  ♦ C0N(ll)^ALPA^RMBaC0N(12)^RMAB  ♦ 
CON ( 13) aRMBaALPAaa^tCON ( 14) aALPAaRMABaCON ( IS) aALPA*8TAaa2* 
C0N(l6)aRMAaBTAaa2  a QaTEMPY  ♦ CON (26) aRaBTAaC0N(27) apagTA 
TEmPZbCON (36) aCON (37) aALPAaCON ( 38 ) aRHAaCON ( 39) aALPAa«2aC0N ( <^0 ) aRHB 

1 aCON (41 ) aBT Aaa2aC0N (42) aALPAaRMAaCON (43) aALPAaRMBaCON (44) aRHAB 
HZOlaCON (28) aCON (29) aBTAaCON (30) aALPAaBTAaCON (31 ) aBTAaRMAaCON ( 32) a 

lBTA*a3aC0N(33)^BTAaALPA^a2aC0N(34)aBTAaRMBaC0N(35)aBTAaRMAB 

2 aRa (EMPZaCON (45) aQaBTAaCON (46) apaxLPA 
THT APaQaocoS (PHI ) -RaosiN (PHI ) 

PSiPa (R*OCOS (PHI ) aOaOSIN (PHI ) ) /DCOS (THTA) 
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MYT»C0N<63»»THTA 

MYT«MYT*C0N(62)*THTAP 

MZTaC0N(64)»PSI 

MZT«MZT*C0N<6l)«PSIP 

MX0I«MX0I*MYT«0SIN<PSI»»0C0S(THTA)*R*1/«Y1-MZT*0SIN<THTA)«RX1/RZ1 
MY01«MY0I*MYT»0C0S (PSl ) •OCOS  «PHI ) ♦ (MYT«0S1N<PSI > *DSIN (THTA) ♦ 

1 MZT»0C0S ( THTA) *RYl/RZl ) •OSIN (PHI ) 
MZOI»MZOI-MYT»DCOS<PSI)»DSIN(PHI)*RZ1/RY1*(MYT*DSIN(PSI)» 

I OSlN (THTA I •RZI/RYl*MZT»DCOS (THTA) ) *OCOS (PHI > 

RETURN 

END 


C 


C 

C 

C 


SUBROUTINE  TEMP(PSItTHTAfPHI,Pl,P2fP3fK) 

IMPLICIT  REAL*8(A-HfO-Z) 

REALMS  HAG 

COMMON/TCOEF/AtB*MA6*ALPA»BTA*OLTA2*RMA»DMA*RMB»DMBtRMABtOMAe 
G AMAsOCOS ( PS I ) •OCOS ( THT A ) 

BTA«-0S1N(PSI>*0C0S<PHI)«0C0S(PSI>*0SIN( THTA) •OSIN (PHI) 

AUPA«0SIN (PSI ) •OSIN (PHI) *0005 (PSD •OSIN (THTA) •OCOS (PHI ) 

0LTA2* (OSIN (PSI ) ) ••Z* (OCOS (PSI ) •OSIN( THTA) ) ••Z 
A»(P3-PI*OSIN(THTA))*OCOS(PSI)^OSIN(PHI)*(PI-P3^DSIN(TMTA)) 

1 ♦OSIN (PSI ) •OCOS (PHI ) ♦PZ^OCOS (PSI ) •OCOS (THTA) •OCOS (PHI ) 

B» (PI^OSIN (THTA) -P3I •OCOS (PSI ) •OCOS (PHI ) ♦ (Pl-P3^0SIN (THTA) ) 

1 ♦OSIN (PSI ) •OSIN (PHI ) ♦PZ^OCOS (PSI ) •OCOS(THTA) ♦OSIN (PHI ) 

MAG-2. 0*0« (PZ^OCOS (PSI ) •OSIN (THTA) ♦P3^0C0S (THTA) •OSIN (PSI ) ) * 

1 OCOS(PSI)^DCOS(THTA) 

G—P3^DS1N  (PSI ) •OCOS  (THTA)  -PZ-DCOS  (PSD  •OSIN  (THTA) 

IF(K.EO.O)  GO  TO  10 

AlaALPA 

BlaBTA 

ALPAaDATAN(Al/6AMA) 

6TA-0ATAN(B1/6AMA) 

Aa ( (OCOS ( ALPA) ) ••Z) * ( A-Al^G/GAHA) /GAMA 
Bb( (DCOS(BTA) )^^2)^(B«B1^6/6AHA)/GAMA 
DLTA2aALPA*^Z  ♦ BTA*^2 
MAG-Z.D^O* (ALPA*A«BTA^B) 

10  CONTINUE 

RMA-DABS(ALPA) 

RMB«DABS(BTA) 

RMABbDABS ( ALPA^BTA) 

OMA>0.0«0 

DMB-0.0^0 

DMABb0.0»0 

1F|ALPA.NE.O.O«0)DMA«ALPA^A/RMA 

IF (ALPA. NE. 0.0*0. AND. BTA.NE,0.D*0)DMAB  • ALPA^BTA^ (BTA*A*ALPA*A) 
1/RMAB 

lF(BTA.NE.0.D«0)OMBaBTA*B/RMB 

RETURN 

END 


SUBROUTINE  FCT(T»X*DERX) 
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IMPLICIT  REAL*8<A-H*0-Z) 

REAL*8  MX01»MY0ItHZ0ItMl*M2,M3»JX«JY«JZtIXP«lYPtlZPfMPM«IPM 
l*lX«lYtIZ 

DIMENSION  X(l)tDERX(I) 

C0MM0N/CaEF/C0N(72) tCCON<72) «NCtNC0N<72) 

COMMON/CONST/JX«JY«JZ*IXP(lYPtIZPtK«IX*IY*IZ 

PHl-XCl) 

THTAbX(2) 

PSI«X(3> 

P»X<A> 

QaX<5) 

MbxI6> 

CAUL  TEHP(PSlfTHTA«PHI*0.O«0f0,D*0»0.O«0(K) 

CALL  MBAR(MXOl»HYOI»MZOl»PSItTHTAtPHl»P»Q*R) 

RJ«l.O*0-JX*JZ 

DERX ( 1 ) bP^OTAN (THTA) • (0»0SIN (PHI ) «R*OCOS (PHI  1 1 
DERX 12) ■Q«OCOS(PHI) -R*DS1N (PHI ) 

DERX (3) ■ (R*DCOS (PHI ) •Q*OSIN (PHI ) ) /OCOS (THTA) 
0ERX(4)s(mX0I«NZ0I»JX-0*R*(JZ*JX«IXP)« 

1 P*Q*(ltO*0-I2P)*JX)/RJ 
DERX (S) -MYOI- (P**2-R»*2) *JY-R*P*I YP 
DERX (6) B (MZ01*MX0I»JZ-P*Q* ( IZP-JX*JZ) • 

1 Q*R*(1.0*0«IXP)*JZ)/RJ 

NE0*(NC*l)*6 
DO  S lB7tNEQ 
S OERX(I)so.O»0 
DO  10  I*1«NC 
IF(NC0N(I).E0*0)G0  TO  10 
P1bX(6«I) 

P2aX'(NC*6«l> 

P3BX(2aNC*6*I) 

01aX(3*NC«6*I) 

U2aX(4*NC«6*I) 

Q3aX(5*NC*6«I) 

CALL  TEMP(PSItTHTA»PHItPlfP2tP3«K) 

CALL  HBAHA(HltH2»H3*PSl«THTA«PHItPtQ*R»PlfP2tP3* 

1 Ql«Q2t03tI) 

DERX (6« I ) aQI«P2« IQ»0SIN (PHI ) *RaOCOS (PHI ) ) / (OCOS (THTA) ) **2 

1 ♦DTAN ( THTA) • (Q2»0SIN (PHI ) ♦P1*0*0C0S (PHI ) 

2 *Q3«0C0S(PHI)-P1»R*0S1N(PH1) ) 

DERX (6«NC«I ) aQ2*0C0S (PHI ) -PiaQ*DSIN (PHI ) -Q3«0SIN (PHI ) 

1 -Pl*R*OCOS(PHl) 

DERX (G*2«NC* I ) a (Q3BDC0S (PHI ) -PlapaOSIN (PHI) ♦02*DSIN (PHI ) 

1 »P1*0«OCOS (PHI > ) /OCOS (THTA) * (R*OCOS (PHI ) 

2 •Q»DSIN(PHI) )*P2*DTAN(THTA)/DCOS(THTA) 

DERX (6«3«NC« I ) > (M1«H3*JX- ( JZ*JX* IXP) • (0*Q3»R«02) 

1 «(1.0*0-lZP)*JXa(P*Q2*Q*Ql) )/RJ 
DERX (6*4«NC* I ) «M2-2.0*0*JY» (P*Q1-R*03) - 
1 IYP*(R*01*P*Q3) 

DERX (6*S*NC« I ) a (M3^MiaJZ- ( IZP-JX«JZ ) • (P*Q2»Q*01 ) 

1 -( 1 .D«0*1XP) aJZ* (Q*03«R»Q2) ) /RJ 
DERX (6«3*NC«I ) "DERX (6*3*NC* I ) • ( (DERX (4) -Q*R) • 

1 2*0«0*CON(60)/(IX»IZ)«HZOI/IX«P*Q«(1.D*0-1ZP)/IX) 

2 aoUMieOtl) /RJ 

DERX(6«4*r4C*I)aDERX(6«4*NC»I)-(P**2-R**2)»DUM(60f I)/IY 
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UERX(6«5*NC«I)B0ERX(6«5*NC«n  (0ERX(6)«P*Q)* 

1 2.U«0*CON<60)/(IX*IZ)«MX01/1Z-Q*R*(1.D«0«IXP)/1Z) 

2 •DUM(60«I)/RJ 

10  continue 

RETURN 

END 


SUBROUTINE  REGOUT (ACONtWXX) 

IMPLICIT  REAL*eiA-H«0-ZI 
REAL*8  MXlfMX2tMYltMY2«MZltMZ2 
DIMENSION  AC0N<72)tWXX(72) 
COMMON/COEF/CON(72) iCC0NI72) *NC«NC0N(72) 
COMMON/CONST2/MXl«MX2tMYltMY2fMZl»MZ2tA 
DO  10  I>ltl6 
AC0N(I}«C0N(1)/MY1 
10  WXXiDBWXXtn/MYl 
DO  12  I«28t35 
ACON(I)aCON(n/MZl 
12  MXX(l)sHXX(l)/MZl 
DO  U I«47»52 
ACON(I)aCON(I>/MXl 

14  WXX(I)aMXX(I)/MXl 
DO  15  l>17t27 
AC0N(I)>C0N(I)/MY2 

15  WXX(1)>MXX(1)/MY2 
DO  16  Ia36»46 
AC0N(I)aC0N(I)/MZ2 

16  WXX(I)aMXX(I)/MZ2 
DO  17  Ia65t72 
HXX(I)aiiilXX(l)/A 

17  ACON(I)aCON(I)/A 
DO  18  Ia53t57 
WXX<n>WXX(l)/MX2 

18  ACON(naCON(n/MX2 
AC0N|58}>C0N(58) 

WXX<58)>WXX(S8> 

ACON(59)>CON<59> 

WXX(59)aHXX(S9> 

4CON(60)*CON<60) 

WXX(60)aMXX<60> 

AC0N(61)aC0N(61)/MZ2 
WXX(61)aMXXI61)/MZ2 
AC0N(63)aC0N(63)/HYl 
ACON ( 64 ) aCON ( 64 ) /MZ 1 
»IXX(63)aWXX(63)/HYl 
MXX(64)aWXX(64)/MZl 
ACON ( 62 ) "CON ( 62 ) /MY2 
MXX(62)aMXX(62)/MY2 
RETURN 
END 

C 

c 

SUBROUTINE  REG (ACON) 
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IMPLICIT  REALMS (A-HtO'Z) 
dimension  AC0N(72) 

REALMS  MXltMX2tMYl*MY2tMZltMZ2 
COMMON/COCF/CON(72) iCCON(72) fNCtNC0N<72l 
C0MH0N/C0NST2/MX1 «HX2tMYl *HY2*HZI *HZ2«  A 
DO  10  I"ltl6 

10  CON(I)>ACON(I)*MYl 
DO  12  Ia28«35 

12  CON(I)bACON(I)*MZ1 
DO  U IbA7«S2 

U C0N(nBAC0N(I)*MXl 
DO  15  I>17«27 

15  CONlI)>ACONai*MY2 
DO  16  Ia36«46 

16  C0N(I)BAC0Nin*MZ2 
DO  17  1b6S«72 

17  C0N(n»AC0N(I)»A 
00  18  lB53t57 

18  C0N(I)«AC0N(I>*MX2 
CON (58) bacon (58) 

CON (59) BACON (59) 

CON(60)bACON(60) 

CON(61)aACON(61)*MZ2 
CON (62) BACON (62) •MYZ 
C0N(63>bAC0N(63)*HY1 
CON ( 64  > B ACON { 64  > BHZ 1 
RETURN 

END 
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Table  3.  Subroutine  Descriptions  for  Both  Asymmetric  and 
Symmetric  Programs 


Name 

Description 

NOISE 

Adds  noise  to  generated  3DOF  motion  or  to  experi- 
mental 3DOF  motion  data.  Input  values  to  the  sub- 
routine specify  the  standard  deviation  and  bias  to 
be  added  to 

m m m 

RANDU 

Random  number  generator. 

GAUSS 

Gaussian  noise  evaluated. 

ICS 

Computes  approximate  boundary  conditions  from  the 
input  3D0F  data  to  be  used  in  place  of  Initial 
guesses . 

PART 

Alters  data  set  to  be  used  in  extraction  technique. 

DRKGS 

Runge-Kutta  fourth-order  integration  scheme  for  NDIM 
simultaneous  differential  equations. 

DUM 

Equivalent  to  the  Kroneker-Delta  function. 

INV 

Provides  solution  to  a set  of  simultaneous  linear 
equations.  Also  used  to  define  inverse  of  a matrix. 

INTGRL 

Calling  subprogram  for  integration  of  parametric, 
kinematic,  and  dynamic  equations.  Evaluates  influ- 
ence coefficient  matrix,  residual  influence  matrix, 
and  total  residual. 

MBARA 

Defines  the  values  of  the  partial  derivatives  of  the 
moment  coefficients  with  respect  to  the  parameters. 

MBAR 

Defines  the  moment  coefficients. 

TEMP 

Defines  variables  required  in  moment  coefficients 
and  derivatives  of  moments  with  respect  to  the 
parameters . 

FCT 

Values  of  derivatives  as  used  in  DRKGS. 

REGOUT 

Conversion  of  parameters  to  output  form. 

REG 

Conversion  of  input  parameters  to  internal  code 
dimensions . 
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Table  4.  Input  Variable  Definitions  for  Symmetric  and 
Asymmetric  Programs 

a.  Namelist  Input  Common  to  Both  Symmetric  and  Asymmetric 
Programs  for  Namelist/NAM/ 


Namelist 

Variable 

Type 

Input 

Values 

Meaning 

Units 

NSWTOl 

Integer 

0,  1,  2 

Switch  Indicating  type  of  analysis; 

0,  normal  data  reduction  by  Chapman- 
Kirk  method  as  outlined  In  Section  3.2; 

- 1,  3D0y  data  generated  from  Initial 
values  of  coefficients;  - 2,  3DOF  data 
generated  and  stored  on  disk  for  later 
analysis , 

NTOT 

Integer 

1,000 

Maximum  array  size  of  Input  and  output 
variables;  time,  6,  '!>. 

NC 

Integer 

30  symmetric 

72  asymmetric 

Maximum  number  of  coefficients  allowed 
to  vary. 

NITER 

Integer 

>1 

Maximum  number  of  Iterations  allowed  to 
reach  a converged  solution. 

DTHIN 

Real 

>0 

Integration  time  step  (DIMIN  must  also 
be  less  than  or  equal  to  the  minimum 
time  step  of  the  Input  data  to  be 
analyzed) . 

BOUND 

1 

Real 

>0 

Convergence  bound  on  Iteration;  convex^ 
gence  Is  assumed  when  difference  In  RMS 
deviation  [Eq.  (26)]  for  two  consecu- 
tive Iterations  is  less  than  or  equal 
to  value  of  BOUND. 

AREA 

Real 

>0 

Reference  aerodynamic  area. 

ft® 

DIA 

Real 

>0 

Reference  aerodynamic  length. 

ft 

IX 

Real 

>0 

Axial  moment  of  inertia. 

slugs-f t® 

lY 

Real 

>0 

Transverse  moment  of  Inertia  in  body- 
fixed  y direction. 

slugs-f t® 

RO 

Real 

>0 

Free^stream  density, 

slugs/ft® 

U 

Real 

>0 

Free-stream  velocity. 

f t/sec 
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Table  4.  Continued 
a.  Continued 


Namelist 

Variable 

Type 

Input 

Values 

Meaning 

Units 

AAOON 

Real 

(aj, 

“No-1,  “Nc) 

Initial  values  of  parameters  which  in- 
clude both  aerodynamic  coefficients 
and  boundary  conditions  on  kinematic 
and  dynamic  equations. “ 

NOON 

Integer 

(0  or  1 for 
each  B]^) 

Array  of  integers  corresponding  to 

AACON,  indicating  which  parameters  are 
allowed  to  change;  - 0,  parameter  is 
constant;  - 1,  parameter  is  varied, 

NSESWH 

Integer 

0 or  1 

Switch  allowing  Gaussian  noise  to  be 
added  to  input  data;  - 1,  noise  added; 

••  0,  no  noise  added. 

SPHI 

Real 

>0 

One  standard  deviation  of  Gaussian 
noise  added  to  the  Euler  angle  4 for 
NSESWH  - 1. 

deg 

SIHA 

Real 

>0 

One  standard  deviation  of  Gaussian 
noise  added  to  the  Euler  angle  6 for 
NSESWH  - 1. 

deg 

SPSI 

Real 

>0 

One  standard  deviation  of  Gaussian 
noise  added  to  the  Euler  angle  w for 
NSESWH  - 1. 

, • 

deg 

AMPHl 

Real 

Constant  bias  in  noise  added  to  Euler 
angle  4 for  NSESWH  - 1, 

deg 

AMTHA 

Real 

Constant  bias  in  noise  added  to  Euler 
angle  9 for  NSESWH  • 1. 

deg 

JQIPSI 

Real 

Constant  bias  in  noise  added  to  Euler 
angle  4 for  NSESWH  - 1. 

deg 

SAMPLE 

Real 

>0 

<1 

Fraction  of  data  sample  to  be  fitted. 

NDELTA 

Integer 

>1 

Only  every  KDELTA  points  of  SAMPLE 
data  are  used  in  fit. 

^See  Tables  4b  and  c for  definitions  of  AACOH  for  symmetric  and  asymmetric  trerslons, 
respectively. 
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Table  4.  Continued 
a.  Concluded 


Namelist 

Variable 

Type 

Input 

Values 

Meaning 

— 

Units 

lUPDAT 

Integer 

0 or  1 

Switch  for  update  of  Initial  values  of 
parameters;  = 0,  Initial  values  not 
updated  for  following  shot;  =>  1, 

Initial  values  updated. 

IMLOOK 

1 

1 

Integer 

0 or  ^ 0 

t 

Switch  defining  conditions  of  Input 
data;  - 0,  Initial  conditions 
evaluated  from  data;  0,  input  values 

of  initial  conditions  used  as  given  In 
AAOOH. 

lAMGLE 

1 

Integer 

0 or  j*  0 

Switch  giving  meaning  of  a and  fi;  = 

0,  a and  S In  terms  of  velocity  ratios; 

^ 0,  a and  S assumed  In  terms  of  angle 
of  attack  and  sideslip. 

DTGEM 

Real 

>0 

Tine  step  of  generated  3DOF  data. 

sec 

IGFNO 

Integer 

Group  number  of  generated  or  Input 
data. 

IBLKNO 

Integer 

Block  number  of  generated  or  Input 
data. 

IRCGEX 

Integer 

Beginning  record  number  of  stored, 
generated  3DOF  data  for  NSWTCH  - 2. 

IHCSTH 

Integer 

Beginning  record  number  of  stored  data 
for  NSNTCH  “ 0 where  IRCSTR  < IRCGEN  - 
3 or  IRCSTR  > IRCGEN  ~ 3. 

i 

SIGPHI 

Real 

>0 

Relative  uncertainty  in  d. 

deg 

5IG1HA 

Real 

>0 

Relative  uncertainty  in  3, 

deg 

SIGPSI 

Real 

>0 

Relative  uncertainty  in  w. 

deg 

NSTART 

Integer 

>1 

Starting  point  of  Input  to  be  fitted. 

IZ 

Real 

>0 

Transverse  moment  of  inertia  In  body- 
fixed  z direction  (asynmetric  program 
only) , 

slugs-f t^ 

IZX 

Real 

>0 

Product  of  Inertia  for  asymmetric 
body  with  mirror  symmetry  about  the 
body-fixed  xz  plane  (asymmetric 
program  only) . 

slugs-ft^  ' 
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Table  4.  Continued 

b.  Namelist  Input  for  AACON  in  Symmetric  Program 


Parameter 

Number 

Meaning 

Units 

Hamel  i st 
Variable 

Meaning 

Units 

1 

2 

3“’ 

aCu/ 3ft 

1/rad 

1/rad 

1/rad^ 

16 

17 

18 

Extra 

Extra 

Extra 

4 

V^(2r 

1/rad^ 

19" 

1/rad 

5 

1/rad^ 

20" 

vlfclr 

1/rad 

6 

1/rad^ 

2l" 

V^=|t 

1/rad 

7 

8 

®y®(|D7) 

l/rad^ 

1/rad^ 

22" 

23" 

3Cn/3v|T 

1/rad 

deg 

9 

ar  /a 

®VM2U  ] 

1/rad^ 

24" 

deg 

10'* 

25 

i(0) 

deg 

11 

ar  /a 1 

®\2U^/ 

1/rad^ 

26 

^(0) 

deg/sec 

12 

v(^) 

1/rad^ 

27 

e{0) 

deg 

13 

\ */ 

1/rad 

28 

6(0) 

deg/sec 

14 

15 

Extra 

1/rad^ 

29 

30 

il<(0) 

^■(0) 

deg 

deg/sec 

I 2 2 

^As  used  here,  4 is  defined  as  s = ya  + 3 . 

‘^These  coefficients  define  the  effects  of  tunnel  flow  nonunifonnities 
on  the  body  ir.otion  and  are  discussed  in  Section  5.3. 

‘^This  coefficient  is  a roll  torque  produced  by  the  gas  bearing  and 
is  discussed  in  Section  5.2. 
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Table  4.  Continued 

c.  Namelist  Input  for  AACON  in  Asymmetric  Program 


Paranieter 

Number 

Meaning 

Units 

Nanelist 

Variable 

Meaning 

Units 

1 

^irio 

19 

v»(#) 

1/rad^ 

2 

ay  3a 

1/rad 

20 

>v>(V) 

l/rad-^ 

2 

3ya|al 

1/rad 

21 

1/rad^ 

4 

3C^/aa2 

1/rad^ 

22 

=v>(#) 

1/rad^ 

5 

ay3a|a| 

1/rad^ 

23 

1/rad^ 

6 

aCn/3a^ 

1/rad^ 

24 

1/rad^ 

7 

ayaa^Ul 

1/rad^ 

25 

1/rad^ 

3 

ay3|a| 

1/rad 

26 

1/rad^ 

9 

3Cn,/362 

1/rad^ 

27 

3C  /a/M^ 

®V^(2U  j 

1/rad^ 

10 

1/rad^ 

28 

''no 

11 

3C^/3a|3l 

1/rad^ 

29 

3C^/3B 

1/rad 

^ *1 

A(- 

3C^/5|ael 

1/rad^ 

30 

oC|^/3aB 

1/rad^ 

13 

3C^/3a2|5| 

1/rad^ 

31 

aySBlal 

1/rad^ 

14 

3C^/5a|a;| 

l/rad”* 

32 

3Cn/3S" 

1/rad^ 

15 

a'yaaB^ 

1/rad^ 

33 

2 

3Cj^/3a^e 

1/rad^ 

1C 

ayas^lal 

1/rad^ 

34 

3C^/30|8l 

1/rad^ 

17 

1/rad^ 

35 

3C^/3B|ae| 

1/rad^ 

18 

^V®(2D^) 

' a»' 

l/rad^ 

36 

’'""ft) 

1/rad 
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Table  4.  Concluded 
c.  Concluded 


Units 

1/rad^ 

1/rad^ 

1/rad^ 

slugs-ft 

1/rad 

1/rad 

1/rad 

1/rad 

deg 

deg 

deg 

deg 

deg 

deg/ sec 
deg/sec 

deg/ sec 


i 


2 
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Table  5. 


Number 
of  Points 
per  Cycle 


Effect  of  Number  of  Cycles  and  Number  of  Points  per 
Cycle  on  C^  q RMS  Deviation  in  Percent  of  C^  q 


Number  of  Cycles  < N 


3Va 

7 

10 

15 

7 

67.0 

21.2 

10.9 

1.0 

15 

46.0 

14.6 

8.2 

5.4 

30 

31.3 

10.9 

5.9 

3.8 

60 

22,0 

7.8 

4.9 

0.5 
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Table  6.  Effect  of  Measurement  Noise  on  Linear  and  Nonlinear 
Coefficients  and  Their  RMS  Deviations 


u« 

deg 

^•ma» 

1/rad 

( Cna)  , 
1/rad 

» 

1/rad 

ECCffltj) , 
1/rad 

mil  6 

l/rad-* 

*'iiici6?5  I 

l/rnd^ 

c ^ 

mq6  ! 
1/rad^ 

0.01 

-0.2360 

0.47  (lO"''"’) 

-4.001 

0.0017 

5.K10H 

0 . 0008 

98.67 

0.43 

0.02 

-0.2360 

0.95  ao~®) 

-3 .998 

0.0033 

5.8108 

0 . 001  6 

97.44 

0.86 

o.o:» 

-0.2360 

0.14  (10-'*) 

-4.004 

0 . 0050 

5.809 

0.0024 

99.59 

1.29 

0.04 

-o,2;if>o 

0.18  (10-4) 

-.3.994 

0 . 0063 

5.814 

0.0031 

97  .32 

1.63 

0.05 

-0.2360 

0.24  (iir'*) 

-4.007 

0.0083 

5.809 

0.0041 

99 . 53 

2.15 

Niite:  E(  ) denotes  HHS  deviulicin. 


^ Table  7.  Tabular  Results  of  Sharp  Cone  Analysis  Using  the  Asymmetric 

Program  for  Fixed  Model  Aerodynamics 


Tunnel-Fixed  Procc.ss  Noi.se  fEq.  (38)J 

I 


Damping 

1 

Stif  fnc.sK 

r ^ 

L 

t4. 

l/rad 

Q 

ti. 

1 /rad 

f-3* 

l/rad 

-T* 

deg 

B 

Sum  of 
Residuals 
[Eq.  (32)], 
deg 

Du  ta 

Sampi ed 
(»oo 

FIk-  17), 
sec 

Value 

Deviation 

-0..57  (10“®) 
0.4  (lOr") 

2.24 

0,11 

0.0163 

0.00016 

-0.0056 

0.00015 

-0.134 

0.011 

0 . 1 52 
O.Ol.'i 

mm 

4.40 

Value 

Dcviiilion 

-0.10  (lor^) 
0.1  (10-6) 

0.0018 

0 . 0003 

1 .35 
. 0.10 

-0.45 

O.OH 

0.0154 

0 . 0001 1 

-0.0054 

0.00014 

-0.118 

0.013 

0.159 

0.014 

0.1158 

8.44 
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Table  8.  Tabular  Results  of  Sphere  Cone  Analysis  Using  the  Asymmetric 
Program  for  Fixed  Values  of  Cgp  = -0.0027  1/rad,  ti  = 0.017 
1/rad,  t4  = 0,  0T  = 0,  i/zj  “ 0 


^ma’ 

1 a 1 ' 

^maipr 

^npla  r 

^'raq' 

1/rad 

l/rad^ 

l/rad^ 

1/rad 

1/rad^ 

1/rad^ 

1/rad 

Value 

-0.240 

-1.87 

-1 .22 

0.265 

0.927 

1.57 

-1.55 

Deviation 

0.005 

0.11 

0.16 

0.005 

0,17 

0.09 

0.64 

Value 

-0.189 

-2.70 

-2.02 

0.269 

0.801 

1.53 

1.25 

Deviation 

0.005 

0.13 

0.18 

0.004 

0.18 

0.10 

0.58 

Value 

-0.172 

-3.06 

-2.20 

0.264 

0.728 

1.62 

1.25 

Deviation 

0.004 

0.12 

0.19 

0.002 

0.19 

0.10 

0.50 

c , 

nr' 

1/rad 

^i!o 

X 10® 

^zx’ 

2 

slugs>-f  t 

H 

4' 

l/rad 

Sum  of 
Residuals 
[Eq,  (32)], 
deg 

Data 

Sampled 

(see 

Fig.  18), 
sec 

Value 

Deviation 

-4.54 

0.63 

3.5 

0.39 

0.0027 

0.0006 

-3.65 

0.20 

-0.0016 

0 . 0004 

0.0775 

5.16 

Value 

Deviation 

-12.3 

0,60 

-1.2 

0.2 

0.0032 

0.0006 

0.62 

0.20 

0.0026 

0.0005 

0.130 

8.32 

Value 

Deviation 

-14.47 

0.57 

-3.5 

0.09 

0.0028 

0 . 0005 

2.72 

0.17 

0.0032 

0.0005 

0.174 

11.98 
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APPENDIX  A 

MATHEMATICAL  DEVELOPMENT:  ASYMMETRIC  VERSION 

The  parameter  identification  technique  as  developed  in  Section  3.0  considered  the 
3DOF  motion  of  axially  symmetric  bodies.  In  order  to  provide  the  greatest  versatility 
in  the  improved  analysis  technique  as  developed  in  this  report,  an  analogous  version 
considering  asymmetric  bodies  has  been  formulated.  The  asymmetric  analysis  considers 
the  3 DOF  motion  of  bodies  with  mirror  plane  symmetry  about  the  x,  z plane  (see  Fig. 
1).  The  governing  equations  as  derived  from  Eqs.  (1)  and  (2)  in  the  body-fixed  coordinate 
system  are 

Jf  = - pq(i;  - J^.IJ  - qrd  i;) 

q = itij.  - (,,2  _ ,.2^  _ ,.p|' 

.Ip  = I - qr(J^.l.,  + i;)  I-  pq(l  - I'J 

where  the  following  definitions  have  been  used 


I ; - « •, 

- I,)/\ 

III 

X 

'< 

U - ('x 

- 'J/I, 

- KJh 

>> 

III 

- IJ/I. 

(A-2) 

with  Mx  = Mx/lx,  My  = My/Iy,  and  M,.  = Mz/I,.. 

The  kinematic  relations  relating  the  Euler  angles  to  the  body-fixed  angular  rates  are  given 
in  Eq.  (A-3). 


0 = p - lull  d (q  sin  <;&  - r cos 

(5  = q cos  <f)  - \-  sin  ^ (A-3) 

lii  -=  {r  cos  {f>  \ q sin  <^)/cos  0 

Integration  of  Eqs.  (A-1)  and  (A-3)  yields  the  3DOF  angular  motion  solution  for  an 
asymmetric  body  assuming  that  the  moment  expansions  as  functions  of  the  Euler  angles 
and  angular  rates  are  given.  For  mirror  symmetry  about  the  x,  z plane  in  the  body-fixed 
axis  system,  the  moment  coefficient  expansions  in  terms  of  a and  j3  as  given  in  Eqs. 
(11)  and  (14)  arc  utilized.  The  variables  a and  arc,  however,  re-evaluated  in  the  body-fixed 
axis  system  leading  to  Eq.  (A-4). 
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P = ''  = -sin  cos  + cos  w sin  0 sin 

^ C»C 

(A-4) 

a = 1!L  = sin  (if  sin  (f>  cos  t!/  sin  0 cos  ti 

r 

30 

The  parameter  extraction  method  as  discussed  in  Section  3.2  is  applied  in  an  analogous 
manner  to  the  asymmetric  analysis.  The  parametric  differential  equations  required  to 
evaluate  the  influence  coefficients  as  used  in  Eq.  (21)  are  derived  from  Eqs.  (A-1)  and 
(A-3)  and  are  shown  in  Eq.  (A-5). 

Pli  = Qj;  + Pgi  sec^  0 (q  sin  o + r cos  6) 

f tan  0 (02i  sin  i + Pjj  q cos  cos  ^ 

-P  ] j r sin  4)) 


hi  ■■= 

Qgi  cos 

6 - P]i 

q sin  9 

- Qgj  sin  6 

- P 1 1 r cos 

hi  = 

(Qjj  cos 

cf>  - P,j 

r sin  ^ 

-(  Qai  sin  4> 

+ Pj  j q cos  4>)  sec  0 

+ (r  cos 

. <j>  + q sin  (^)  Pai  tan  0 see  0 

JQ,i 

= Mn  - 

Msplx  - 

i;)  (qQ3i 

- rQa;!  (A-5) 

t- 

(1  - r,).i. 

(p  Q2i  + 

Mai  - J 

y(2pQ,i 

- 2rQ3 

,i)  - i;(rQ, 

i ^ pQa.) 

i%i 

= -''3i 

M,i.I,  - 

(i;  - J 

J,)  (pQai 

•-  qQji) 

- 

(qQa.  - rQa;) 

where  Ph,  P2i,  Psi,  Mn,  and  M2i,  and  M3i  are  defined  as  before  and  Q]  i = 3p/9a,, 
Q2i  = 9q/9ai,  Qsi  = 3r/0ai.  The  initial  conditions  for  the  parametric,  kinematic,  and 
dynamic  equations  as  given  in  Eqs.  (A-1),  (A-3),  and  (A-5)  are  listed  below. 


o(0)  = a.,^ 

p (0)  - a.yj  _ 2 

c 

c 

0(0)  = £1^  __  1 

q(0)  = a„  _ , 

‘ c 

c 

0(0)  = _3 

r(0)  = 

c 

(A-6) 


174 


AEDC-TR-78-10 


C 

PgjCO)  = Sx|^_4  j 

Q2i(0)  = j 

(A-7) 

PjitO)  = 5|^,  _3  j 
c 

QgilO)  = 

Integration  of  the  parametric  and  angular  motion  equations  and  evaluation  of  the 
corrections  to  the  parameters,  Uj,  are  determined  in  a manner  similar  to  that  discussed 
in  Section  3.2  for  the  axially  symmetric  analysis. 

In  addition  to  the  moment  coefficients  as  given  in  Eq.  (14)  and  the  boundary 
conditions,  the  product  of  inertia,  has  also  been  included  in  the  list  of  parameters, 
ai,  to  be  determined.  This  option  vi'as  provided  in  order  to  circumvent  the  lack  of  simple 
techniques  for  determining  the  product  of  inertia  for  asymmetric  configurations.  In  Ref. 
21,  values  of  I^^  as  large  as  6 percent  of  ly  have  been  extracted  successfully  from  generated 
motion  with  measurement  noise. 
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^lo  > 

Bi 

1 bill-" 


Cij 


CR 

Cfio 

Cfip 

Cfiq 

Qr 

Cfii 

C)2(j 
Cm 
Cm  p 
Cmq 

Cm  q S 2 

Cmr 

Cm  a 
Cm  a 

Cm  a I o I 


NOMENCLATURE 

Aerodynamic  coefficients  to  be  extracted,  see  Section  3.2 

Taylor  series  expansion  coefficients  as  used  in  Eqs.  (12),  (14),  and 
(16)  (i  = 1,  m,  n) 

Diagonal  components  of  the  inverse  matrix  of  Qj. 

Taylor  series  expansion  coefficients  as  used  in  Eqs.  (12),  (14),  and 
(16)  (i=  l,m,'n) 

Nf  X Nc  matrix  of  influence  coefficients,  Eq.  (24) 
Rolling-moment  coefficient,  Eq.  (7) 

Trim  roll  torque  coefficient,  see  Eq.  (16) 

9Cg/a(dp/2UJ,  1/rad 
aCg/9(dq/2UJ,  1/rad 
9Cg/9(dr/2UJ,  1/rad 
9Cg/9(da/2UJ,  1/rad 
9Ce/9(dP/2UJ,  1/rad 
Pitching-moment  coefficient,  Eq.  (7) 

9Cm/9(dp/2UJ,  1/rad 
9Cm/3(dq/2UJ,  1/rad 
9Cm/3(dq52/2UJ,  l/rad3 
9Cm/3(dr/2UJ,  1/rad 
9Cm/3as  l/tad 
9C^/9(da/2UJ,  1/rad 
9C„/9alal,  l/rad^ 
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Cm  a I 0 i 

Cm  aS  ^ 
Cm  p 

C„ 


Cnp 

Cnq 

Cfi  r 

Cna 

Cjip 

Cnp’ 

Cnp  I a I 
Cnp  Ip  I 


Cio> 

Dj 

d 

d3>  di,  ds 

dN 

E 


eo 

ei,  62 

fios 


3Cm/3al|3l,  l/rad^ 

9Cm/3a52,  l/rad^ 

3Cm/3(d/3/2UJ.  1/rad 
Yawing-moment  coefficient,  Eq.  (7) 

3Cn/3(dp/2UJ,  1/rad 
3Cn/3(dq/2UJ,  1/rad 
3Cn/3(dr/2UJ,  1/rad 
3C„/3(da/2UJ,  1/rad 

3Cn/3/3,  1/rad 
3C„/3(d/3/2UJ,  1/rad 

3C„/3plal,  l/rad2 

3Cn/3|3l|3l,  l/rad2 

Taylor  series  expansion  coefficients  as  used  in  Eqs.  (12),  (14), 
and  (16)  (i  = 1,  m,  n) 

Residual  influence  matrix,  Eq.  (25) 

Reference  length  (base  diameter),  ft 

Coefficients  as  used  in  Eq.  (16) 

Nose  diameter,  ft 

RMS  deviation  between  experimental  data  and  converged  solution, 
rad 

Trim  roll  coefficient,  see  Eq.  (16) 

Proportionality  constants  relating  a,  q and  /?,  r 

Taylor  series  expansion  coefficients  as  used  in  Eqs.  (12),  (14), 
and  (16)  (i  = 1,  m,  n) 
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S3  > S4  > Ss 


I 

lx 

Ixx  O*"  ^x 
Ixy 

ly 

ly  y or  ly 

lyz 

i; 

Izx 

Izz 

^ t 

1,  J,  k 


J 

Jx 

Jy 

Jz 

fi 


Mii,  M2i,  Msi 
k^x  > > ^z 


Coefficients  as  used  in  Eq.  (16) 

Moment  of  momentum  vector  defined  in  Eq.  (2) 

Components  of  moment  of  momentum  vector 

Transverse  moment  of  inertia  for  axially  symmetric  body,  slugs-ft^ 

Dimensionless  moment  of  inertia,  Eq.  (A-2) 

Moment  of  inertia  about  x axis,  Fig.  1,  slugs-ft^ 

Product  of  inertia  in  x,  y plane,  Fig.  1,  slugs-ft^ 

Dimensionless  moment  of  inertia,  Eq.  (A-2) 

Moment  of  inertia  about  y axis,  Fig.  1,  slugs-ft^ 

Product  of  inertia  in  y,  z plane,  Fig.  1,  slugs-ft^ 

Dimensionless  moment  of  inertia,  Eq.  (A-2) 

Product  of  inertia  in  x,  z plane.  Fig.  1,  slugs-ft^ 

Moment  of  inertia  about  z axis,  Fig.  1,  slugs-ft^ 

Unit  vectors  along  coordinate  axis,  Fig.  1 
1 - JxJz 

Dimensionless  product  of  inertia,  Eq.  (A-2) 

Dimensionless  product  of  inertia,  Eq.  (A-2) 

Dimensionless  product  of  inertia,  Eq.  (A-2) 

Distance  between  bearing  pivot  point  and  mass  center,  see  Fig. 
13,  ft 

Moment  vector  acting  at  mass  eg  [see  Eq.  (1)  and  discussion] 
Equal  to  aMj^/Bai,  aMy/dai,  and  bUjdzi,  respectively 
Components  of  moment  vector 
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MyT 

MzT 

N 

No 

n 

Pli>  P2i.  Psi 


Pli2,  P2iK,  P3i£ 

P 

Qli.  Q2i.  Qsi 
q 

q- 

R 

Ri«,  R22,  Rsfi 
r 
S 
t 


tl,  t2,  t3,  U 

tm  ax 

u 

u£ 


Equal  to  Mx/lx.  My /I,  and  Mz/I,  respectively,  for  the  symmetric 
program  and  Mx/Ix,  My/Iy,  and  respectively,  for  the 

asymmetric  program,  Ibf/slug-ft 

Tunnel-fixed  moment,  Eq.  (38),  ft-lbf 

Tunnel-fixed  moment,  Eq.  (38),  ft-lbf 

Total  number  of  time  points  in  fit 

Total  number  of  coefficients  to  be  extracted 

Number  of  iterations  per  case 

Equal  to  30/9ai,  dfl/Saj,  and  9\(//9ai,  respectively 

Value  of  Pii,  P2i,  and  Psj,  respectively,  at  time  point  fi 

Roll  velocity,  deg/sec 

Equal  to  3p/3ai,  3q/3ai,  3r/3ai,  respectively 
Pitch  velocity,  deg/sec 
Dynamic  pressure  of  free-stream,  psfa 
Dimensionless  inertia  ratio,  Ix/I 

Residuals  in  0,  0,  and  i//,  respectively,  at  each  time  point  £,  deg 

Yaw  rate,  deg/sec 

Reference  area  (base  area),  ft^ 

Time,  sec 

Tunnel-fixed  moment  coefficients  as  defined  in  Eq.  (38),  1/rad 

Data  sample  time,  sec 

Free-stream  velocity,  ft/sec 

Velocity  along  x axis,  see  Fig.  1,  ft/sec 

Residual  in  roll,  deg 
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V 

VK 

W 

w 

WJJ 

X 

X 

x', 

Y 

y 

z 

z 

a 

Oa 

Ov 


Aai 

At 

5 

h 


Velocity  along  y axis,  see  Fig.  1,  ft/sec 

Residual  in  yaw,  deg 

Pendulum  weight,  Ibf 

Velocity  along  z axis,  see  Fig.  1,  ft/sec 

Residual  in  pitch,  deg 

Tunnel-fixed  axis,  see  Fig.  2 

Body-fixed  axis,  see  Fig.  1 

Rotated  axis  system  as  defined  in  Fig.  3 

Tunnel-fixed  axis,  see  Fig.  2 

Body-fixed  axis,  see  Fig.  1 

Tunnel-fixed  axis,  see  Fig.  2 

Body-fixed  axis,  see  Fig.  1 

Velocity  ratio  as  given  in  Eq.  (9)  or  Eq.  (A-4) 

An^e  of  attack  as  defined  in  Eq.  (33),  rad 

Velocity  ratio  as  used  by  Tobak,  Eq.  (32) 

Velocity  ratio  as  gjven  in  Eq.  (9)  or  Eq.  (A-4) 

Angle  of  attack  as  defined  in  Eq.  (33),  rad 

Velocity  ratio  as  used  by  Tobak,  Eq.  (32) 

Corrections  to  the  coefficients  ai 

Integration  time  interval,  sec 

Equal  to  y/a^  + rad 

Kroneker  Delta  function 
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5ai 


e 

6c 

Oj 


6i 

69.0 

a 

0(f>,  0$,  0^ 


<l> 

09 
09o 
0m  9 
0 
0T 

09 

09o 


0m  9 


£2 


RMS  deviation  in  extracted  coefficient  aj 
Pitch  angle,  see  Fig.  2 
Cone  semiangle 

Artificial  tunnel-fixed  flow  angularity  correction,  pitch,  deg 

Calculated  pitch  angle  at  time  point  )2,  deg 

Value  of  69  for  the  given  constants,  ai,  deg 

Measured  pitch  angle  at  time  point  £,  deg 

One  standard  deviation,  Gaussian  noise,  deg 

Standard  deviation  of  relative  uncertainty  in  (f>,  0,  and  0 
measurement,  respectively,  deg 

Computer  run  time  of  codes,  Eq.  (31),  sec 

Roll  angle,  deg 

Calculated  roll  angle  at  time  point  K,  deg 
Value  of  09  for  the  given  constants,  aj,  deg 
Measured  roll  angle  at  time  point  £,  deg 
Yaw  angle,  see  Fig.  2,  deg 

Artificial  tunnel-fixed  flow  angularity  correction,  yaw,  deg 

Calculated  yaw  angle  at  time  point  C,  deg 

Value  of  09  for  the  given  constants,  aj,  deg 

Measured  yaw  angle  at  time  point  fi,  deg 

Sting  roll  position  as  used  in  Eq.  (36),  deg 
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Superscripts 

^ Coefficients  in  plane-fixed  or  aeroballistic  axis  system 

d/dt 


Subscript 

o 


Initial  value  of  coefficient 


182 


